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