qdivrem.c revision 92889
11573Srgrimes/*- 21573Srgrimes * Copyright (c) 1992, 1993 31573Srgrimes * The Regents of the University of California. All rights reserved. 41573Srgrimes * 51573Srgrimes * This software was developed by the Computer Systems Engineering group 61573Srgrimes * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and 71573Srgrimes * contributed to Berkeley. 81573Srgrimes * 91573Srgrimes * Redistribution and use in source and binary forms, with or without 101573Srgrimes * modification, are permitted provided that the following conditions 111573Srgrimes * are met: 121573Srgrimes * 1. Redistributions of source code must retain the above copyright 131573Srgrimes * notice, this list of conditions and the following disclaimer. 141573Srgrimes * 2. Redistributions in binary form must reproduce the above copyright 151573Srgrimes * notice, this list of conditions and the following disclaimer in the 161573Srgrimes * documentation and/or other materials provided with the distribution. 171573Srgrimes * 3. All advertising materials mentioning features or use of this software 181573Srgrimes * must display the following acknowledgement: 191573Srgrimes * This product includes software developed by the University of 201573Srgrimes * California, Berkeley and its contributors. 211573Srgrimes * 4. Neither the name of the University nor the names of its contributors 221573Srgrimes * may be used to endorse or promote products derived from this software 231573Srgrimes * without specific prior written permission. 241573Srgrimes * 251573Srgrimes * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 261573Srgrimes * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 271573Srgrimes * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 281573Srgrimes * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 291573Srgrimes * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 301573Srgrimes * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 311573Srgrimes * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 321573Srgrimes * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 331573Srgrimes * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 341573Srgrimes * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 351573Srgrimes * SUCH DAMAGE. 361573Srgrimes */ 371573Srgrimes 381573Srgrimes#if defined(LIBC_SCCS) && !defined(lint) 391573Srgrimesstatic char sccsid[] = "@(#)qdivrem.c 8.1 (Berkeley) 6/4/93"; 401573Srgrimes#endif /* LIBC_SCCS and not lint */ 4192889Sobrien#include <sys/cdefs.h> 4292889Sobrien__FBSDID("$FreeBSD: head/lib/libc/quad/qdivrem.c 92889 2002-03-21 18:49:23Z obrien $"); 431573Srgrimes 441573Srgrimes/* 451573Srgrimes * Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed), 461573Srgrimes * section 4.3.1, pp. 257--259. 471573Srgrimes */ 481573Srgrimes 491573Srgrimes#include "quad.h" 501573Srgrimes 511573Srgrimes#define B (1 << HALF_BITS) /* digit base */ 521573Srgrimes 531573Srgrimes/* Combine two `digits' to make a single two-digit number. */ 541573Srgrimes#define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b)) 551573Srgrimes 561573Srgrimes/* select a type for digits in base B: use unsigned short if they fit */ 571573Srgrimes#if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff 581573Srgrimestypedef unsigned short digit; 591573Srgrimes#else 601573Srgrimestypedef u_long digit; 611573Srgrimes#endif 621573Srgrimes 631573Srgrimes/* 641573Srgrimes * Shift p[0]..p[len] left `sh' bits, ignoring any bits that 651573Srgrimes * `fall out' the left (there never will be any such anyway). 661573Srgrimes * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS. 671573Srgrimes */ 681573Srgrimesstatic void 6992889Sobrienshl(digit *p, int len, int sh) 701573Srgrimes{ 7192889Sobrien int i; 721573Srgrimes 731573Srgrimes for (i = 0; i < len; i++) 741573Srgrimes p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh)); 751573Srgrimes p[i] = LHALF(p[i] << sh); 761573Srgrimes} 771573Srgrimes 781573Srgrimes/* 791573Srgrimes * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v. 801573Srgrimes * 811573Srgrimes * We do this in base 2-sup-HALF_BITS, so that all intermediate products 821573Srgrimes * fit within u_long. As a consequence, the maximum length dividend and 831573Srgrimes * divisor are 4 `digits' in this base (they are shorter if they have 841573Srgrimes * leading zeros). 851573Srgrimes */ 861573Srgrimesu_quad_t 871573Srgrimes__qdivrem(uq, vq, arq) 881573Srgrimes u_quad_t uq, vq, *arq; 891573Srgrimes{ 901573Srgrimes union uu tmp; 911573Srgrimes digit *u, *v, *q; 9292889Sobrien digit v1, v2; 931573Srgrimes u_long qhat, rhat, t; 941573Srgrimes int m, n, d, j, i; 951573Srgrimes digit uspace[5], vspace[5], qspace[5]; 961573Srgrimes 971573Srgrimes /* 981573Srgrimes * Take care of special cases: divide by zero, and u < v. 991573Srgrimes */ 1001573Srgrimes if (vq == 0) { 1011573Srgrimes /* divide by zero. */ 1021573Srgrimes static volatile const unsigned int zero = 0; 1031573Srgrimes 1041573Srgrimes tmp.ul[H] = tmp.ul[L] = 1 / zero; 1051573Srgrimes if (arq) 1061573Srgrimes *arq = uq; 1071573Srgrimes return (tmp.q); 1081573Srgrimes } 1091573Srgrimes if (uq < vq) { 1101573Srgrimes if (arq) 1111573Srgrimes *arq = uq; 1121573Srgrimes return (0); 1131573Srgrimes } 1141573Srgrimes u = &uspace[0]; 1151573Srgrimes v = &vspace[0]; 1161573Srgrimes q = &qspace[0]; 1171573Srgrimes 1181573Srgrimes /* 1191573Srgrimes * Break dividend and divisor into digits in base B, then 1201573Srgrimes * count leading zeros to determine m and n. When done, we 1211573Srgrimes * will have: 1221573Srgrimes * u = (u[1]u[2]...u[m+n]) sub B 1231573Srgrimes * v = (v[1]v[2]...v[n]) sub B 1241573Srgrimes * v[1] != 0 1251573Srgrimes * 1 < n <= 4 (if n = 1, we use a different division algorithm) 1261573Srgrimes * m >= 0 (otherwise u < v, which we already checked) 1271573Srgrimes * m + n = 4 1281573Srgrimes * and thus 1291573Srgrimes * m = 4 - n <= 2 1301573Srgrimes */ 1311573Srgrimes tmp.uq = uq; 1321573Srgrimes u[0] = 0; 1331573Srgrimes u[1] = HHALF(tmp.ul[H]); 1341573Srgrimes u[2] = LHALF(tmp.ul[H]); 1351573Srgrimes u[3] = HHALF(tmp.ul[L]); 1361573Srgrimes u[4] = LHALF(tmp.ul[L]); 1371573Srgrimes tmp.uq = vq; 1381573Srgrimes v[1] = HHALF(tmp.ul[H]); 1391573Srgrimes v[2] = LHALF(tmp.ul[H]); 1401573Srgrimes v[3] = HHALF(tmp.ul[L]); 1411573Srgrimes v[4] = LHALF(tmp.ul[L]); 1421573Srgrimes for (n = 4; v[1] == 0; v++) { 1431573Srgrimes if (--n == 1) { 1441573Srgrimes u_long rbj; /* r*B+u[j] (not root boy jim) */ 1451573Srgrimes digit q1, q2, q3, q4; 1461573Srgrimes 1471573Srgrimes /* 1481573Srgrimes * Change of plan, per exercise 16. 1491573Srgrimes * r = 0; 1501573Srgrimes * for j = 1..4: 1511573Srgrimes * q[j] = floor((r*B + u[j]) / v), 1521573Srgrimes * r = (r*B + u[j]) % v; 1531573Srgrimes * We unroll this completely here. 1541573Srgrimes */ 1551573Srgrimes t = v[2]; /* nonzero, by definition */ 1561573Srgrimes q1 = u[1] / t; 1571573Srgrimes rbj = COMBINE(u[1] % t, u[2]); 1581573Srgrimes q2 = rbj / t; 1591573Srgrimes rbj = COMBINE(rbj % t, u[3]); 1601573Srgrimes q3 = rbj / t; 1611573Srgrimes rbj = COMBINE(rbj % t, u[4]); 1621573Srgrimes q4 = rbj / t; 1631573Srgrimes if (arq) 1641573Srgrimes *arq = rbj % t; 1651573Srgrimes tmp.ul[H] = COMBINE(q1, q2); 1661573Srgrimes tmp.ul[L] = COMBINE(q3, q4); 1671573Srgrimes return (tmp.q); 1681573Srgrimes } 1691573Srgrimes } 1701573Srgrimes 1711573Srgrimes /* 1721573Srgrimes * By adjusting q once we determine m, we can guarantee that 1731573Srgrimes * there is a complete four-digit quotient at &qspace[1] when 1741573Srgrimes * we finally stop. 1751573Srgrimes */ 1761573Srgrimes for (m = 4 - n; u[1] == 0; u++) 1771573Srgrimes m--; 1781573Srgrimes for (i = 4 - m; --i >= 0;) 1791573Srgrimes q[i] = 0; 1801573Srgrimes q += 4 - m; 1811573Srgrimes 1821573Srgrimes /* 1831573Srgrimes * Here we run Program D, translated from MIX to C and acquiring 1841573Srgrimes * a few minor changes. 1851573Srgrimes * 1861573Srgrimes * D1: choose multiplier 1 << d to ensure v[1] >= B/2. 1871573Srgrimes */ 1881573Srgrimes d = 0; 1891573Srgrimes for (t = v[1]; t < B / 2; t <<= 1) 1901573Srgrimes d++; 1911573Srgrimes if (d > 0) { 1921573Srgrimes shl(&u[0], m + n, d); /* u <<= d */ 1931573Srgrimes shl(&v[1], n - 1, d); /* v <<= d */ 1941573Srgrimes } 1951573Srgrimes /* 1961573Srgrimes * D2: j = 0. 1971573Srgrimes */ 1981573Srgrimes j = 0; 1991573Srgrimes v1 = v[1]; /* for D3 -- note that v[1..n] are constant */ 2001573Srgrimes v2 = v[2]; /* for D3 */ 2011573Srgrimes do { 20292889Sobrien digit uj0, uj1, uj2; 2038870Srgrimes 2041573Srgrimes /* 2051573Srgrimes * D3: Calculate qhat (\^q, in TeX notation). 2061573Srgrimes * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and 2071573Srgrimes * let rhat = (u[j]*B + u[j+1]) mod v[1]. 2081573Srgrimes * While rhat < B and v[2]*qhat > rhat*B+u[j+2], 2091573Srgrimes * decrement qhat and increase rhat correspondingly. 2101573Srgrimes * Note that if rhat >= B, v[2]*qhat < rhat*B. 2111573Srgrimes */ 2121573Srgrimes uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */ 2131573Srgrimes uj1 = u[j + 1]; /* for D3 only */ 2141573Srgrimes uj2 = u[j + 2]; /* for D3 only */ 2151573Srgrimes if (uj0 == v1) { 2161573Srgrimes qhat = B; 2171573Srgrimes rhat = uj1; 2181573Srgrimes goto qhat_too_big; 2191573Srgrimes } else { 2201573Srgrimes u_long n = COMBINE(uj0, uj1); 2211573Srgrimes qhat = n / v1; 2221573Srgrimes rhat = n % v1; 2231573Srgrimes } 2241573Srgrimes while (v2 * qhat > COMBINE(rhat, uj2)) { 2251573Srgrimes qhat_too_big: 2261573Srgrimes qhat--; 2271573Srgrimes if ((rhat += v1) >= B) 2281573Srgrimes break; 2291573Srgrimes } 2301573Srgrimes /* 2311573Srgrimes * D4: Multiply and subtract. 2321573Srgrimes * The variable `t' holds any borrows across the loop. 2331573Srgrimes * We split this up so that we do not require v[0] = 0, 2341573Srgrimes * and to eliminate a final special case. 2351573Srgrimes */ 2361573Srgrimes for (t = 0, i = n; i > 0; i--) { 2371573Srgrimes t = u[i + j] - v[i] * qhat - t; 2381573Srgrimes u[i + j] = LHALF(t); 2391573Srgrimes t = (B - HHALF(t)) & (B - 1); 2401573Srgrimes } 2411573Srgrimes t = u[j] - t; 2421573Srgrimes u[j] = LHALF(t); 2431573Srgrimes /* 2441573Srgrimes * D5: test remainder. 2451573Srgrimes * There is a borrow if and only if HHALF(t) is nonzero; 2461573Srgrimes * in that (rare) case, qhat was too large (by exactly 1). 2471573Srgrimes * Fix it by adding v[1..n] to u[j..j+n]. 2481573Srgrimes */ 2491573Srgrimes if (HHALF(t)) { 2501573Srgrimes qhat--; 2511573Srgrimes for (t = 0, i = n; i > 0; i--) { /* D6: add back. */ 2521573Srgrimes t += u[i + j] + v[i]; 2531573Srgrimes u[i + j] = LHALF(t); 2541573Srgrimes t = HHALF(t); 2551573Srgrimes } 2561573Srgrimes u[j] = LHALF(u[j] + t); 2571573Srgrimes } 2581573Srgrimes q[j] = qhat; 2591573Srgrimes } while (++j <= m); /* D7: loop on j. */ 2601573Srgrimes 2611573Srgrimes /* 2621573Srgrimes * If caller wants the remainder, we have to calculate it as 2631573Srgrimes * u[m..m+n] >> d (this is at most n digits and thus fits in 2641573Srgrimes * u[m+1..m+n], but we may need more source digits). 2651573Srgrimes */ 2661573Srgrimes if (arq) { 2671573Srgrimes if (d) { 2681573Srgrimes for (i = m + n; i > m; --i) 2691573Srgrimes u[i] = (u[i] >> d) | 2701573Srgrimes LHALF(u[i - 1] << (HALF_BITS - d)); 2711573Srgrimes u[i] = 0; 2721573Srgrimes } 2731573Srgrimes tmp.ul[H] = COMBINE(uspace[1], uspace[2]); 2741573Srgrimes tmp.ul[L] = COMBINE(uspace[3], uspace[4]); 2751573Srgrimes *arq = tmp.q; 2761573Srgrimes } 2771573Srgrimes 2781573Srgrimes tmp.ul[H] = COMBINE(qspace[1], qspace[2]); 2791573Srgrimes tmp.ul[L] = COMBINE(qspace[3], qspace[4]); 2801573Srgrimes return (tmp.q); 2811573Srgrimes} 282