Annotation of sys/lib/libkern/qdivrem.c, Revision 1.1.1.1
1.1 nbrk 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: * 3. 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: #if defined(LIBC_SCCS) && !defined(lint)
35: static char rcsid[] = "$OpenBSD: qdivrem.c,v 1.8 2005/02/13 03:37:14 jsg Exp $";
36: #endif /* LIBC_SCCS and not lint */
37:
38: /*
39: * Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed),
40: * section 4.3.1, pp. 257--259.
41: */
42:
43: #include "quad.h"
44:
45: #define B ((int)1 << HALF_BITS) /* digit base */
46:
47: /* Combine two `digits' to make a single two-digit number. */
48: #define COMBINE(a, b) (((u_int)(a) << HALF_BITS) | (b))
49:
50: /* select a type for digits in base B: use unsigned short if they fit */
51: #if UINT_MAX == 0xffffffffU && USHRT_MAX >= 0xffff
52: typedef unsigned short digit;
53: #else
54: typedef u_int digit;
55: #endif
56:
57: static void shl(digit *p, int len, int sh);
58:
59: /*
60: * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
61: *
62: * We do this in base 2-sup-HALF_BITS, so that all intermediate products
63: * fit within u_int. As a consequence, the maximum length dividend and
64: * divisor are 4 `digits' in this base (they are shorter if they have
65: * leading zeros).
66: */
67: u_quad_t
68: __qdivrem(u_quad_t uq, u_quad_t vq, u_quad_t *arq)
69: {
70: union uu tmp;
71: digit *u, *v, *q;
72: digit v1, v2;
73: u_int qhat, rhat, t;
74: int m, n, d, j, i;
75: digit uspace[5], vspace[5], qspace[5];
76:
77: /*
78: * Take care of special cases: divide by zero, and u < v.
79: */
80: if (vq == 0) {
81: /* divide by zero. */
82: static volatile const unsigned int zero = 0;
83:
84: tmp.ul[H] = tmp.ul[L] = 1 / zero;
85: if (arq)
86: *arq = uq;
87: return (tmp.q);
88: }
89: if (uq < vq) {
90: if (arq)
91: *arq = uq;
92: return (0);
93: }
94: u = &uspace[0];
95: v = &vspace[0];
96: q = &qspace[0];
97:
98: /*
99: * Break dividend and divisor into digits in base B, then
100: * count leading zeros to determine m and n. When done, we
101: * will have:
102: * u = (u[1]u[2]...u[m+n]) sub B
103: * v = (v[1]v[2]...v[n]) sub B
104: * v[1] != 0
105: * 1 < n <= 4 (if n = 1, we use a different division algorithm)
106: * m >= 0 (otherwise u < v, which we already checked)
107: * m + n = 4
108: * and thus
109: * m = 4 - n <= 2
110: */
111: tmp.uq = uq;
112: u[0] = 0;
113: u[1] = (digit)HHALF(tmp.ul[H]);
114: u[2] = (digit)LHALF(tmp.ul[H]);
115: u[3] = (digit)HHALF(tmp.ul[L]);
116: u[4] = (digit)LHALF(tmp.ul[L]);
117: tmp.uq = vq;
118: v[1] = (digit)HHALF(tmp.ul[H]);
119: v[2] = (digit)LHALF(tmp.ul[H]);
120: v[3] = (digit)HHALF(tmp.ul[L]);
121: v[4] = (digit)LHALF(tmp.ul[L]);
122: for (n = 4; v[1] == 0; v++) {
123: if (--n == 1) {
124: u_int rbj; /* r*B+u[j] (not root boy jim) */
125: digit q1, q2, q3, q4;
126:
127: /*
128: * Change of plan, per exercise 16.
129: * r = 0;
130: * for j = 1..4:
131: * q[j] = floor((r*B + u[j]) / v),
132: * r = (r*B + u[j]) % v;
133: * We unroll this completely here.
134: */
135: t = v[2]; /* nonzero, by definition */
136: q1 = (digit)(u[1] / t);
137: rbj = COMBINE(u[1] % t, u[2]);
138: q2 = (digit)(rbj / t);
139: rbj = COMBINE(rbj % t, u[3]);
140: q3 = (digit)(rbj / t);
141: rbj = COMBINE(rbj % t, u[4]);
142: q4 = (digit)(rbj / t);
143: if (arq)
144: *arq = rbj % t;
145: tmp.ul[H] = COMBINE(q1, q2);
146: tmp.ul[L] = COMBINE(q3, q4);
147: return (tmp.q);
148: }
149: }
150:
151: /*
152: * By adjusting q once we determine m, we can guarantee that
153: * there is a complete four-digit quotient at &qspace[1] when
154: * we finally stop.
155: */
156: for (m = 4 - n; u[1] == 0; u++)
157: m--;
158: for (i = 4 - m; --i >= 0;)
159: q[i] = 0;
160: q += 4 - m;
161:
162: /*
163: * Here we run Program D, translated from MIX to C and acquiring
164: * a few minor changes.
165: *
166: * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
167: */
168: d = 0;
169: for (t = v[1]; t < B / 2; t <<= 1)
170: d++;
171: if (d > 0) {
172: shl(&u[0], m + n, d); /* u <<= d */
173: shl(&v[1], n - 1, d); /* v <<= d */
174: }
175: /*
176: * D2: j = 0.
177: */
178: j = 0;
179: v1 = v[1]; /* for D3 -- note that v[1..n] are constant */
180: v2 = v[2]; /* for D3 */
181: do {
182: digit uj0, uj1, uj2;
183:
184: /*
185: * D3: Calculate qhat (\^q, in TeX notation).
186: * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
187: * let rhat = (u[j]*B + u[j+1]) mod v[1].
188: * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
189: * decrement qhat and increase rhat correspondingly.
190: * Note that if rhat >= B, v[2]*qhat < rhat*B.
191: */
192: uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */
193: uj1 = u[j + 1]; /* for D3 only */
194: uj2 = u[j + 2]; /* for D3 only */
195: if (uj0 == v1) {
196: qhat = B;
197: rhat = uj1;
198: goto qhat_too_big;
199: } else {
200: u_int nn = COMBINE(uj0, uj1);
201: qhat = nn / v1;
202: rhat = nn % v1;
203: }
204: while (v2 * qhat > COMBINE(rhat, uj2)) {
205: qhat_too_big:
206: qhat--;
207: if ((rhat += v1) >= B)
208: break;
209: }
210: /*
211: * D4: Multiply and subtract.
212: * The variable `t' holds any borrows across the loop.
213: * We split this up so that we do not require v[0] = 0,
214: * and to eliminate a final special case.
215: */
216: for (t = 0, i = n; i > 0; i--) {
217: t = u[i + j] - v[i] * qhat - t;
218: u[i + j] = (digit)LHALF(t);
219: t = (B - HHALF(t)) & (B - 1);
220: }
221: t = u[j] - t;
222: u[j] = (digit)LHALF(t);
223: /*
224: * D5: test remainder.
225: * There is a borrow if and only if HHALF(t) is nonzero;
226: * in that (rare) case, qhat was too large (by exactly 1).
227: * Fix it by adding v[1..n] to u[j..j+n].
228: */
229: if (HHALF(t)) {
230: qhat--;
231: for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
232: t += u[i + j] + v[i];
233: u[i + j] = (digit)LHALF(t);
234: t = HHALF(t);
235: }
236: u[j] = (digit)LHALF(u[j] + t);
237: }
238: q[j] = (digit)qhat;
239: } while (++j <= m); /* D7: loop on j. */
240:
241: /*
242: * If caller wants the remainder, we have to calculate it as
243: * u[m..m+n] >> d (this is at most n digits and thus fits in
244: * u[m+1..m+n], but we may need more source digits).
245: */
246: if (arq) {
247: if (d) {
248: for (i = m + n; i > m; --i)
249: u[i] = (digit)(((u_int)u[i] >> d) |
250: LHALF((u_int)u[i - 1] << (HALF_BITS - d)));
251: u[i] = 0;
252: }
253: tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
254: tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
255: *arq = tmp.q;
256: }
257:
258: tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
259: tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
260: return (tmp.q);
261: }
262:
263: /*
264: * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
265: * `fall out' the left (there never will be any such anyway).
266: * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS.
267: */
268: static void
269: shl(digit *p, int len, int sh)
270: {
271: int i;
272:
273: for (i = 0; i < len; i++)
274: p[i] = (digit)(LHALF((u_int)p[i] << sh) |
275: ((u_int)p[i + 1] >> (HALF_BITS - sh)));
276: p[i] = (digit)(LHALF((u_int)p[i] << sh));
277: }
CVSweb
![[BACK]](/icons/back.gif){kind=icon}
Return to qdivrem.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / lib / libkern