xref: /illumos-gate/usr/src/boot/libsa/qdivrem.c (revision 22028508)
1199767f8SToomas Soome /*-
2199767f8SToomas Soome  * Copyright (c) 1992, 1993
3199767f8SToomas Soome  *	The Regents of the University of California.  All rights reserved.
4199767f8SToomas Soome  *
5199767f8SToomas Soome  * This software was developed by the Computer Systems Engineering group
6199767f8SToomas Soome  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
7199767f8SToomas Soome  * contributed to Berkeley.
8199767f8SToomas Soome  *
9199767f8SToomas Soome  * Redistribution and use in source and binary forms, with or without
10199767f8SToomas Soome  * modification, are permitted provided that the following conditions
11199767f8SToomas Soome  * are met:
12199767f8SToomas Soome  * 1. Redistributions of source code must retain the above copyright
13199767f8SToomas Soome  *    notice, this list of conditions and the following disclaimer.
14199767f8SToomas Soome  * 2. Redistributions in binary form must reproduce the above copyright
15199767f8SToomas Soome  *    notice, this list of conditions and the following disclaimer in the
16199767f8SToomas Soome  *    documentation and/or other materials provided with the distribution.
171974da4bSToomas Soome  * 3. Neither the name of the University nor the names of its contributors
18199767f8SToomas Soome  *    may be used to endorse or promote products derived from this software
19199767f8SToomas Soome  *    without specific prior written permission.
20199767f8SToomas Soome  *
21199767f8SToomas Soome  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22199767f8SToomas Soome  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23199767f8SToomas Soome  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24199767f8SToomas Soome  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25199767f8SToomas Soome  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26199767f8SToomas Soome  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27199767f8SToomas Soome  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28199767f8SToomas Soome  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29199767f8SToomas Soome  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30199767f8SToomas Soome  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31199767f8SToomas Soome  * SUCH DAMAGE.
32199767f8SToomas Soome  *
33199767f8SToomas Soome  * 	From: Id: qdivrem.c,v 1.7 1997/11/07 09:20:40 phk Exp
34199767f8SToomas Soome  */
35199767f8SToomas Soome 
36199767f8SToomas Soome #include <sys/cdefs.h>
37d8b4f927SToomas Soome #include <sys/stddef.h>
38199767f8SToomas Soome 
39199767f8SToomas Soome /*
40199767f8SToomas Soome  * Multiprecision divide.  This algorithm is from Knuth vol. 2 (2nd ed),
41199767f8SToomas Soome  * section 4.3.1, pp. 257--259.
42199767f8SToomas Soome  */
43199767f8SToomas Soome 
44199767f8SToomas Soome #include "quad.h"
45199767f8SToomas Soome 
46199767f8SToomas Soome #define	B	(1 << HALF_BITS)	/* digit base */
47199767f8SToomas Soome 
48199767f8SToomas Soome /* Combine two `digits' to make a single two-digit number. */
49199767f8SToomas Soome #define	COMBINE(a, b) (((u_int)(a) << HALF_BITS) | (b))
50199767f8SToomas Soome 
51199767f8SToomas Soome _Static_assert(sizeof(int) / 2 == sizeof(short),
52199767f8SToomas Soome 	"Bitwise functions in libstand are broken on this architecture\n");
53199767f8SToomas Soome 
54199767f8SToomas Soome /* select a type for digits in base B: use unsigned short if they fit */
55199767f8SToomas Soome typedef unsigned short digit;
56199767f8SToomas Soome 
57199767f8SToomas Soome /*
58199767f8SToomas Soome  * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
59199767f8SToomas Soome  * `fall out' the left (there never will be any such anyway).
60199767f8SToomas Soome  * We may assume len >= 0.  NOTE THAT THIS WRITES len+1 DIGITS.
61199767f8SToomas Soome  */
62199767f8SToomas Soome static void
shl(digit * p,int len,int sh)63199767f8SToomas Soome shl(digit *p, int len, int sh)
64199767f8SToomas Soome {
65199767f8SToomas Soome 	int i;
66199767f8SToomas Soome 
67199767f8SToomas Soome 	for (i = 0; i < len; i++)
68199767f8SToomas Soome 		p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh));
69199767f8SToomas Soome 	p[i] = LHALF(p[i] << sh);
70199767f8SToomas Soome }
71199767f8SToomas Soome 
72199767f8SToomas Soome /*
73d8b4f927SToomas Soome  * __udivmoddi4(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
74199767f8SToomas Soome  *
75199767f8SToomas Soome  * We do this in base 2-sup-HALF_BITS, so that all intermediate products
76199767f8SToomas Soome  * fit within u_int.  As a consequence, the maximum length dividend and
77199767f8SToomas Soome  * divisor are 4 `digits' in this base (they are shorter if they have
78199767f8SToomas Soome  * leading zeros).
79199767f8SToomas Soome  */
80199767f8SToomas Soome u_quad_t
__udivmoddi4(u_quad_t uq,u_quad_t vq,u_quad_t * arq)81d8b4f927SToomas Soome __udivmoddi4(u_quad_t uq, u_quad_t vq, u_quad_t *arq)
82199767f8SToomas Soome {
83199767f8SToomas Soome 	union uu tmp;
84199767f8SToomas Soome 	digit *u, *v, *q;
85199767f8SToomas Soome 	digit v1, v2;
86199767f8SToomas Soome 	u_int qhat, rhat, t;
87199767f8SToomas Soome 	int m, n, d, j, i;
88199767f8SToomas Soome 	digit uspace[5], vspace[5], qspace[5];
89199767f8SToomas Soome 
90199767f8SToomas Soome 	/*
91199767f8SToomas Soome 	 * Take care of special cases: divide by zero, and u < v.
92199767f8SToomas Soome 	 */
93199767f8SToomas Soome 	if (vq == 0) {
94199767f8SToomas Soome 		/* divide by zero. */
95199767f8SToomas Soome 		static volatile const unsigned int zero = 0;
96199767f8SToomas Soome 
97199767f8SToomas Soome 		tmp.ul[H] = tmp.ul[L] = 1 / zero;
98199767f8SToomas Soome 		if (arq)
99199767f8SToomas Soome 			*arq = uq;
100199767f8SToomas Soome 		return (tmp.q);
101199767f8SToomas Soome 	}
102199767f8SToomas Soome 	if (uq < vq) {
103199767f8SToomas Soome 		if (arq)
104199767f8SToomas Soome 			*arq = uq;
105199767f8SToomas Soome 		return (0);
106199767f8SToomas Soome 	}
107199767f8SToomas Soome 	u = &uspace[0];
108199767f8SToomas Soome 	v = &vspace[0];
109199767f8SToomas Soome 	q = &qspace[0];
110199767f8SToomas Soome 
111199767f8SToomas Soome 	/*
112199767f8SToomas Soome 	 * Break dividend and divisor into digits in base B, then
113199767f8SToomas Soome 	 * count leading zeros to determine m and n.  When done, we
114199767f8SToomas Soome 	 * will have:
115199767f8SToomas Soome 	 *	u = (u[1]u[2]...u[m+n]) sub B
116199767f8SToomas Soome 	 *	v = (v[1]v[2]...v[n]) sub B
117199767f8SToomas Soome 	 *	v[1] != 0
118199767f8SToomas Soome 	 *	1 < n <= 4 (if n = 1, we use a different division algorithm)
119199767f8SToomas Soome 	 *	m >= 0 (otherwise u < v, which we already checked)
120199767f8SToomas Soome 	 *	m + n = 4
121199767f8SToomas Soome 	 * and thus
122199767f8SToomas Soome 	 *	m = 4 - n <= 2
123199767f8SToomas Soome 	 */
124199767f8SToomas Soome 	tmp.uq = uq;
125199767f8SToomas Soome 	u[0] = 0;
126199767f8SToomas Soome 	u[1] = HHALF(tmp.ul[H]);
127199767f8SToomas Soome 	u[2] = LHALF(tmp.ul[H]);
128199767f8SToomas Soome 	u[3] = HHALF(tmp.ul[L]);
129199767f8SToomas Soome 	u[4] = LHALF(tmp.ul[L]);
130199767f8SToomas Soome 	tmp.uq = vq;
131199767f8SToomas Soome 	v[1] = HHALF(tmp.ul[H]);
132199767f8SToomas Soome 	v[2] = LHALF(tmp.ul[H]);
133199767f8SToomas Soome 	v[3] = HHALF(tmp.ul[L]);
134199767f8SToomas Soome 	v[4] = LHALF(tmp.ul[L]);
135199767f8SToomas Soome 	for (n = 4; v[1] == 0; v++) {
136199767f8SToomas Soome 		if (--n == 1) {
137199767f8SToomas Soome 			u_int rbj;	/* r*B+u[j] (not root boy jim) */
138199767f8SToomas Soome 			digit q1, q2, q3, q4;
139199767f8SToomas Soome 
140199767f8SToomas Soome 			/*
141199767f8SToomas Soome 			 * Change of plan, per exercise 16.
142199767f8SToomas Soome 			 *	r = 0;
143199767f8SToomas Soome 			 *	for j = 1..4:
144199767f8SToomas Soome 			 *		q[j] = floor((r*B + u[j]) / v),
145199767f8SToomas Soome 			 *		r = (r*B + u[j]) % v;
146199767f8SToomas Soome 			 * We unroll this completely here.
147199767f8SToomas Soome 			 */
148199767f8SToomas Soome 			t = v[2];	/* nonzero, by definition */
149199767f8SToomas Soome 			q1 = u[1] / t;
150199767f8SToomas Soome 			rbj = COMBINE(u[1] % t, u[2]);
151199767f8SToomas Soome 			q2 = rbj / t;
152199767f8SToomas Soome 			rbj = COMBINE(rbj % t, u[3]);
153199767f8SToomas Soome 			q3 = rbj / t;
154199767f8SToomas Soome 			rbj = COMBINE(rbj % t, u[4]);
155199767f8SToomas Soome 			q4 = rbj / t;
156199767f8SToomas Soome 			if (arq)
157199767f8SToomas Soome 				*arq = rbj % t;
158199767f8SToomas Soome 			tmp.ul[H] = COMBINE(q1, q2);
159199767f8SToomas Soome 			tmp.ul[L] = COMBINE(q3, q4);
160199767f8SToomas Soome 			return (tmp.q);
161199767f8SToomas Soome 		}
162199767f8SToomas Soome 	}
163199767f8SToomas Soome 
164199767f8SToomas Soome 	/*
165199767f8SToomas Soome 	 * By adjusting q once we determine m, we can guarantee that
166199767f8SToomas Soome 	 * there is a complete four-digit quotient at &qspace[1] when
167199767f8SToomas Soome 	 * we finally stop.
168199767f8SToomas Soome 	 */
169199767f8SToomas Soome 	for (m = 4 - n; u[1] == 0; u++)
170199767f8SToomas Soome 		m--;
171199767f8SToomas Soome 	for (i = 4 - m; --i >= 0;)
172199767f8SToomas Soome 		q[i] = 0;
173199767f8SToomas Soome 	q += 4 - m;
174199767f8SToomas Soome 
175199767f8SToomas Soome 	/*
176199767f8SToomas Soome 	 * Here we run Program D, translated from MIX to C and acquiring
177199767f8SToomas Soome 	 * a few minor changes.
178199767f8SToomas Soome 	 *
179199767f8SToomas Soome 	 * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
180199767f8SToomas Soome 	 */
181199767f8SToomas Soome 	d = 0;
182199767f8SToomas Soome 	for (t = v[1]; t < B / 2; t <<= 1)
183199767f8SToomas Soome 		d++;
184199767f8SToomas Soome 	if (d > 0) {
185199767f8SToomas Soome 		shl(&u[0], m + n, d);		/* u <<= d */
186199767f8SToomas Soome 		shl(&v[1], n - 1, d);		/* v <<= d */
187199767f8SToomas Soome 	}
188199767f8SToomas Soome 	/*
189199767f8SToomas Soome 	 * D2: j = 0.
190199767f8SToomas Soome 	 */
191199767f8SToomas Soome 	j = 0;
192199767f8SToomas Soome 	v1 = v[1];	/* for D3 -- note that v[1..n] are constant */
193199767f8SToomas Soome 	v2 = v[2];	/* for D3 */
194199767f8SToomas Soome 	do {
195199767f8SToomas Soome 		digit uj0, uj1, uj2;
196199767f8SToomas Soome 
197199767f8SToomas Soome 		/*
198199767f8SToomas Soome 		 * D3: Calculate qhat (\^q, in TeX notation).
199199767f8SToomas Soome 		 * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
200199767f8SToomas Soome 		 * let rhat = (u[j]*B + u[j+1]) mod v[1].
201199767f8SToomas Soome 		 * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
202199767f8SToomas Soome 		 * decrement qhat and increase rhat correspondingly.
203199767f8SToomas Soome 		 * Note that if rhat >= B, v[2]*qhat < rhat*B.
204199767f8SToomas Soome 		 */
205199767f8SToomas Soome 		uj0 = u[j + 0];	/* for D3 only -- note that u[j+...] change */
206199767f8SToomas Soome 		uj1 = u[j + 1];	/* for D3 only */
207199767f8SToomas Soome 		uj2 = u[j + 2];	/* for D3 only */
208199767f8SToomas Soome 		if (uj0 == v1) {
209199767f8SToomas Soome 			qhat = B;
210199767f8SToomas Soome 			rhat = uj1;
211199767f8SToomas Soome 			goto qhat_too_big;
212199767f8SToomas Soome 		} else {
213199767f8SToomas Soome 			u_int nn = COMBINE(uj0, uj1);
214199767f8SToomas Soome 			qhat = nn / v1;
215199767f8SToomas Soome 			rhat = nn % v1;
216199767f8SToomas Soome 		}
217199767f8SToomas Soome 		while (v2 * qhat > COMBINE(rhat, uj2)) {
218199767f8SToomas Soome 	qhat_too_big:
219199767f8SToomas Soome 			qhat--;
220199767f8SToomas Soome 			if ((rhat += v1) >= B)
221199767f8SToomas Soome 				break;
222199767f8SToomas Soome 		}
223199767f8SToomas Soome 		/*
224199767f8SToomas Soome 		 * D4: Multiply and subtract.
225199767f8SToomas Soome 		 * The variable `t' holds any borrows across the loop.
226199767f8SToomas Soome 		 * We split this up so that we do not require v[0] = 0,
227199767f8SToomas Soome 		 * and to eliminate a final special case.
228199767f8SToomas Soome 		 */
229199767f8SToomas Soome 		for (t = 0, i = n; i > 0; i--) {
230199767f8SToomas Soome 			t = u[i + j] - v[i] * qhat - t;
231199767f8SToomas Soome 			u[i + j] = LHALF(t);
232199767f8SToomas Soome 			t = (B - HHALF(t)) & (B - 1);
233199767f8SToomas Soome 		}
234199767f8SToomas Soome 		t = u[j] - t;
235199767f8SToomas Soome 		u[j] = LHALF(t);
236199767f8SToomas Soome 		/*
237199767f8SToomas Soome 		 * D5: test remainder.
238199767f8SToomas Soome 		 * There is a borrow if and only if HHALF(t) is nonzero;
239199767f8SToomas Soome 		 * in that (rare) case, qhat was too large (by exactly 1).
240199767f8SToomas Soome 		 * Fix it by adding v[1..n] to u[j..j+n].
241199767f8SToomas Soome 		 */
242199767f8SToomas Soome 		if (HHALF(t)) {
243199767f8SToomas Soome 			qhat--;
244199767f8SToomas Soome 			for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
245199767f8SToomas Soome 				t += u[i + j] + v[i];
246199767f8SToomas Soome 				u[i + j] = LHALF(t);
247199767f8SToomas Soome 				t = HHALF(t);
248199767f8SToomas Soome 			}
249199767f8SToomas Soome 			u[j] = LHALF(u[j] + t);
250199767f8SToomas Soome 		}
251199767f8SToomas Soome 		q[j] = qhat;
252199767f8SToomas Soome 	} while (++j <= m);		/* D7: loop on j. */
253199767f8SToomas Soome 
254199767f8SToomas Soome 	/*
255199767f8SToomas Soome 	 * If caller wants the remainder, we have to calculate it as
256199767f8SToomas Soome 	 * u[m..m+n] >> d (this is at most n digits and thus fits in
257199767f8SToomas Soome 	 * u[m+1..m+n], but we may need more source digits).
258199767f8SToomas Soome 	 */
259199767f8SToomas Soome 	if (arq) {
260199767f8SToomas Soome 		if (d) {
261199767f8SToomas Soome 			for (i = m + n; i > m; --i)
262199767f8SToomas Soome 				u[i] = (u[i] >> d) |
263199767f8SToomas Soome 				    LHALF(u[i - 1] << (HALF_BITS - d));
264199767f8SToomas Soome 			u[i] = 0;
265199767f8SToomas Soome 		}
266199767f8SToomas Soome 		tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
267199767f8SToomas Soome 		tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
268199767f8SToomas Soome 		*arq = tmp.q;
269199767f8SToomas Soome 	}
270199767f8SToomas Soome 
271199767f8SToomas Soome 	tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
272199767f8SToomas Soome 	tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
273199767f8SToomas Soome 	return (tmp.q);
274199767f8SToomas Soome }
275199767f8SToomas Soome 
276199767f8SToomas Soome /*
277199767f8SToomas Soome  * Divide two unsigned quads.
278199767f8SToomas Soome  */
279199767f8SToomas Soome 
280199767f8SToomas Soome u_quad_t
__udivdi3(u_quad_t a,u_quad_t b)281d8b4f927SToomas Soome __udivdi3(u_quad_t a, u_quad_t b)
282199767f8SToomas Soome {
283199767f8SToomas Soome 
284d8b4f927SToomas Soome 	return (__udivmoddi4(a, b, NULL));
285199767f8SToomas Soome }
286199767f8SToomas Soome 
287199767f8SToomas Soome /*
288199767f8SToomas Soome  * Return remainder after dividing two unsigned quads.
289199767f8SToomas Soome  */
290199767f8SToomas Soome u_quad_t
__umoddi3(u_quad_t a,u_quad_t b)291d8b4f927SToomas Soome __umoddi3(u_quad_t a, u_quad_t b)
292199767f8SToomas Soome {
293199767f8SToomas Soome 	u_quad_t r;
294199767f8SToomas Soome 
295d8b4f927SToomas Soome 	(void)__udivmoddi4(a, b, &r);
296199767f8SToomas Soome 	return (r);
297199767f8SToomas Soome }
298199767f8SToomas Soome 
299199767f8SToomas Soome /*
300199767f8SToomas Soome  * Divide two signed quads.
301199767f8SToomas Soome  * ??? if -1/2 should produce -1 on this machine, this code is wrong
302199767f8SToomas Soome  */
303199767f8SToomas Soome quad_t
__divdi3(quad_t a,quad_t b)304d8b4f927SToomas Soome __divdi3(quad_t a, quad_t b)
305199767f8SToomas Soome {
306199767f8SToomas Soome 	u_quad_t ua, ub, uq;
307199767f8SToomas Soome 	int neg;
308199767f8SToomas Soome 
309199767f8SToomas Soome 	if (a < 0)
310199767f8SToomas Soome 		ua = -(u_quad_t)a, neg = 1;
311199767f8SToomas Soome 	else
312199767f8SToomas Soome 		ua = a, neg = 0;
313199767f8SToomas Soome 	if (b < 0)
314199767f8SToomas Soome 		ub = -(u_quad_t)b, neg ^= 1;
315199767f8SToomas Soome 	else
316199767f8SToomas Soome 		ub = b;
317d8b4f927SToomas Soome 	uq = __udivmoddi4(ua, ub, NULL);
318199767f8SToomas Soome 	return (neg ? -uq : uq);
319199767f8SToomas Soome }
320199767f8SToomas Soome 
321199767f8SToomas Soome /*
322199767f8SToomas Soome  * Return remainder after dividing two signed quads.
323199767f8SToomas Soome  *
324199767f8SToomas Soome  * XXX
325199767f8SToomas Soome  * If -1/2 should produce -1 on this machine, this code is wrong.
326199767f8SToomas Soome  */
327199767f8SToomas Soome quad_t
__moddi3(quad_t a,quad_t b)328d8b4f927SToomas Soome __moddi3(quad_t a, quad_t b)
329199767f8SToomas Soome {
330199767f8SToomas Soome 	u_quad_t ua, ub, ur;
331199767f8SToomas Soome 	int neg;
332199767f8SToomas Soome 
333199767f8SToomas Soome 	if (a < 0)
334199767f8SToomas Soome 		ua = -(u_quad_t)a, neg = 1;
335199767f8SToomas Soome 	else
336199767f8SToomas Soome 		ua = a, neg = 0;
337199767f8SToomas Soome 	if (b < 0)
338199767f8SToomas Soome 		ub = -(u_quad_t)b;
339199767f8SToomas Soome 	else
340199767f8SToomas Soome 		ub = b;
341d8b4f927SToomas Soome 	(void)__udivmoddi4(ua, ub, &ur);
342199767f8SToomas Soome 	return (neg ? -ur : ur);
343199767f8SToomas Soome }
344d8b4f927SToomas Soome 
345d8b4f927SToomas Soome quad_t
__divmoddi4(quad_t a,quad_t b,quad_t * r)346d8b4f927SToomas Soome __divmoddi4(quad_t a, quad_t b, quad_t *r)
347d8b4f927SToomas Soome {
348d8b4f927SToomas Soome 	quad_t d;
349d8b4f927SToomas Soome 
350d8b4f927SToomas Soome 	d = __divdi3(a, b);
351402902c3SToomas Soome 	if (r != NULL)
352402902c3SToomas Soome 		*r = a - (b * d);
353d8b4f927SToomas Soome 
354d8b4f927SToomas Soome 	return (d);
355d8b4f927SToomas Soome }
356