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