xref: /freebsd/sys/libkern/qdivrem.c (revision 1e413cf93298b5b97441a21d9a50fdcd0ee9945e)
1 /*-
2  * Copyright (c) 1992, 1993
3  *	The Regents of the University of California.  All rights reserved.
4  *
5  * This software was developed by the Computer Systems Engineering group
6  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
7  * contributed to Berkeley.
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions
11  * are met:
12  * 1. Redistributions of source code must retain the above copyright
13  *    notice, this list of conditions and the following disclaimer.
14  * 2. Redistributions in binary form must reproduce the above copyright
15  *    notice, this list of conditions and the following disclaimer in the
16  *    documentation and/or other materials provided with the distribution.
17  * 4. Neither the name of the University nor the names of its contributors
18  *    may be used to endorse or promote products derived from this software
19  *    without specific prior written permission.
20  *
21  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31  * SUCH DAMAGE.
32  */
33 
34 #include <sys/cdefs.h>
35 __FBSDID("$FreeBSD$");
36 
37 /*
38  * Multiprecision divide.  This algorithm is from Knuth vol. 2 (2nd ed),
39  * section 4.3.1, pp. 257--259.
40  */
41 
42 #include <libkern/quad.h>
43 
44 #define	B	(1 << HALF_BITS)	/* digit base */
45 
46 /* Combine two `digits' to make a single two-digit number. */
47 #define	COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b))
48 
49 /* select a type for digits in base B: use unsigned short if they fit */
50 #if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff
51 typedef unsigned short digit;
52 #else
53 typedef u_long digit;
54 #endif
55 
56 /*
57  * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
58  * `fall out' the left (there never will be any such anyway).
59  * We may assume len >= 0.  NOTE THAT THIS WRITES len+1 DIGITS.
60  */
61 static void
62 shl(register digit *p, register int len, register int sh)
63 {
64 	register int i;
65 
66 	for (i = 0; i < len; i++)
67 		p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh));
68 	p[i] = LHALF(p[i] << sh);
69 }
70 
71 /*
72  * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
73  *
74  * We do this in base 2-sup-HALF_BITS, so that all intermediate products
75  * fit within u_long.  As a consequence, the maximum length dividend and
76  * divisor are 4 `digits' in this base (they are shorter if they have
77  * leading zeros).
78  */
79 u_quad_t
80 __qdivrem(uq, vq, arq)
81 	u_quad_t uq, vq, *arq;
82 {
83 	union uu tmp;
84 	digit *u, *v, *q;
85 	register digit v1, v2;
86 	u_long qhat, rhat, t;
87 	int m, n, d, j, i;
88 	digit uspace[5], vspace[5], qspace[5];
89 
90 	/*
91 	 * Take care of special cases: divide by zero, and u < v.
92 	 */
93 	if (vq == 0) {
94 		/* divide by zero. */
95 		static volatile const unsigned int zero = 0;
96 
97 		tmp.ul[H] = tmp.ul[L] = 1 / zero;
98 		if (arq)
99 			*arq = uq;
100 		return (tmp.q);
101 	}
102 	if (uq < vq) {
103 		if (arq)
104 			*arq = uq;
105 		return (0);
106 	}
107 	u = &uspace[0];
108 	v = &vspace[0];
109 	q = &qspace[0];
110 
111 	/*
112 	 * Break dividend and divisor into digits in base B, then
113 	 * count leading zeros to determine m and n.  When done, we
114 	 * will have:
115 	 *	u = (u[1]u[2]...u[m+n]) sub B
116 	 *	v = (v[1]v[2]...v[n]) sub B
117 	 *	v[1] != 0
118 	 *	1 < n <= 4 (if n = 1, we use a different division algorithm)
119 	 *	m >= 0 (otherwise u < v, which we already checked)
120 	 *	m + n = 4
121 	 * and thus
122 	 *	m = 4 - n <= 2
123 	 */
124 	tmp.uq = uq;
125 	u[0] = 0;
126 	u[1] = HHALF(tmp.ul[H]);
127 	u[2] = LHALF(tmp.ul[H]);
128 	u[3] = HHALF(tmp.ul[L]);
129 	u[4] = LHALF(tmp.ul[L]);
130 	tmp.uq = vq;
131 	v[1] = HHALF(tmp.ul[H]);
132 	v[2] = LHALF(tmp.ul[H]);
133 	v[3] = HHALF(tmp.ul[L]);
134 	v[4] = LHALF(tmp.ul[L]);
135 	for (n = 4; v[1] == 0; v++) {
136 		if (--n == 1) {
137 			u_long rbj;	/* r*B+u[j] (not root boy jim) */
138 			digit q1, q2, q3, q4;
139 
140 			/*
141 			 * Change of plan, per exercise 16.
142 			 *	r = 0;
143 			 *	for j = 1..4:
144 			 *		q[j] = floor((r*B + u[j]) / v),
145 			 *		r = (r*B + u[j]) % v;
146 			 * We unroll this completely here.
147 			 */
148 			t = v[2];	/* nonzero, by definition */
149 			q1 = u[1] / t;
150 			rbj = COMBINE(u[1] % t, u[2]);
151 			q2 = rbj / t;
152 			rbj = COMBINE(rbj % t, u[3]);
153 			q3 = rbj / t;
154 			rbj = COMBINE(rbj % t, u[4]);
155 			q4 = rbj / t;
156 			if (arq)
157 				*arq = rbj % t;
158 			tmp.ul[H] = COMBINE(q1, q2);
159 			tmp.ul[L] = COMBINE(q3, q4);
160 			return (tmp.q);
161 		}
162 	}
163 
164 	/*
165 	 * By adjusting q once we determine m, we can guarantee that
166 	 * there is a complete four-digit quotient at &qspace[1] when
167 	 * we finally stop.
168 	 */
169 	for (m = 4 - n; u[1] == 0; u++)
170 		m--;
171 	for (i = 4 - m; --i >= 0;)
172 		q[i] = 0;
173 	q += 4 - m;
174 
175 	/*
176 	 * Here we run Program D, translated from MIX to C and acquiring
177 	 * a few minor changes.
178 	 *
179 	 * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
180 	 */
181 	d = 0;
182 	for (t = v[1]; t < B / 2; t <<= 1)
183 		d++;
184 	if (d > 0) {
185 		shl(&u[0], m + n, d);		/* u <<= d */
186 		shl(&v[1], n - 1, d);		/* v <<= d */
187 	}
188 	/*
189 	 * D2: j = 0.
190 	 */
191 	j = 0;
192 	v1 = v[1];	/* for D3 -- note that v[1..n] are constant */
193 	v2 = v[2];	/* for D3 */
194 	do {
195 		register digit uj0, uj1, uj2;
196 
197 		/*
198 		 * D3: Calculate qhat (\^q, in TeX notation).
199 		 * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
200 		 * let rhat = (u[j]*B + u[j+1]) mod v[1].
201 		 * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
202 		 * decrement qhat and increase rhat correspondingly.
203 		 * Note that if rhat >= B, v[2]*qhat < rhat*B.
204 		 */
205 		uj0 = u[j + 0];	/* for D3 only -- note that u[j+...] change */
206 		uj1 = u[j + 1];	/* for D3 only */
207 		uj2 = u[j + 2];	/* for D3 only */
208 		if (uj0 == v1) {
209 			qhat = B;
210 			rhat = uj1;
211 			goto qhat_too_big;
212 		} else {
213 			u_long nn = COMBINE(uj0, uj1);
214 			qhat = nn / v1;
215 			rhat = nn % v1;
216 		}
217 		while (v2 * qhat > COMBINE(rhat, uj2)) {
218 	qhat_too_big:
219 			qhat--;
220 			if ((rhat += v1) >= B)
221 				break;
222 		}
223 		/*
224 		 * D4: Multiply and subtract.
225 		 * The variable `t' holds any borrows across the loop.
226 		 * We split this up so that we do not require v[0] = 0,
227 		 * and to eliminate a final special case.
228 		 */
229 		for (t = 0, i = n; i > 0; i--) {
230 			t = u[i + j] - v[i] * qhat - t;
231 			u[i + j] = LHALF(t);
232 			t = (B - HHALF(t)) & (B - 1);
233 		}
234 		t = u[j] - t;
235 		u[j] = LHALF(t);
236 		/*
237 		 * D5: test remainder.
238 		 * There is a borrow if and only if HHALF(t) is nonzero;
239 		 * in that (rare) case, qhat was too large (by exactly 1).
240 		 * Fix it by adding v[1..n] to u[j..j+n].
241 		 */
242 		if (HHALF(t)) {
243 			qhat--;
244 			for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
245 				t += u[i + j] + v[i];
246 				u[i + j] = LHALF(t);
247 				t = HHALF(t);
248 			}
249 			u[j] = LHALF(u[j] + t);
250 		}
251 		q[j] = qhat;
252 	} while (++j <= m);		/* D7: loop on j. */
253 
254 	/*
255 	 * If caller wants the remainder, we have to calculate it as
256 	 * u[m..m+n] >> d (this is at most n digits and thus fits in
257 	 * u[m+1..m+n], but we may need more source digits).
258 	 */
259 	if (arq) {
260 		if (d) {
261 			for (i = m + n; i > m; --i)
262 				u[i] = (u[i] >> d) |
263 				    LHALF(u[i - 1] << (HALF_BITS - d));
264 			u[i] = 0;
265 		}
266 		tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
267 		tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
268 		*arq = tmp.q;
269 	}
270 
271 	tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
272 	tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
273 	return (tmp.q);
274 }
275