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