qdivrem.c revision 1573
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 */ 411573Srgrimes 421573Srgrimes/* 431573Srgrimes * Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed), 441573Srgrimes * section 4.3.1, pp. 257--259. 451573Srgrimes */ 461573Srgrimes 471573Srgrimes#include "quad.h" 481573Srgrimes 491573Srgrimes#define B (1 << HALF_BITS) /* digit base */ 501573Srgrimes 511573Srgrimes/* Combine two `digits' to make a single two-digit number. */ 521573Srgrimes#define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b)) 531573Srgrimes 541573Srgrimes/* select a type for digits in base B: use unsigned short if they fit */ 551573Srgrimes#if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff 561573Srgrimestypedef unsigned short digit; 571573Srgrimes#else 581573Srgrimestypedef u_long digit; 591573Srgrimes#endif 601573Srgrimes 611573Srgrimes/* 621573Srgrimes * Shift p[0]..p[len] left `sh' bits, ignoring any bits that 631573Srgrimes * `fall out' the left (there never will be any such anyway). 641573Srgrimes * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS. 651573Srgrimes */ 661573Srgrimesstatic void 671573Srgrimesshl(register digit *p, register int len, register int sh) 681573Srgrimes{ 691573Srgrimes register int i; 701573Srgrimes 711573Srgrimes for (i = 0; i < len; i++) 721573Srgrimes p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh)); 731573Srgrimes p[i] = LHALF(p[i] << sh); 741573Srgrimes} 751573Srgrimes 761573Srgrimes/* 771573Srgrimes * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v. 781573Srgrimes * 791573Srgrimes * We do this in base 2-sup-HALF_BITS, so that all intermediate products 801573Srgrimes * fit within u_long. As a consequence, the maximum length dividend and 811573Srgrimes * divisor are 4 `digits' in this base (they are shorter if they have 821573Srgrimes * leading zeros). 831573Srgrimes */ 841573Srgrimesu_quad_t 851573Srgrimes__qdivrem(uq, vq, arq) 861573Srgrimes u_quad_t uq, vq, *arq; 871573Srgrimes{ 881573Srgrimes union uu tmp; 891573Srgrimes digit *u, *v, *q; 901573Srgrimes register digit v1, v2; 911573Srgrimes u_long qhat, rhat, t; 921573Srgrimes int m, n, d, j, i; 931573Srgrimes digit uspace[5], vspace[5], qspace[5]; 941573Srgrimes 951573Srgrimes /* 961573Srgrimes * Take care of special cases: divide by zero, and u < v. 971573Srgrimes */ 981573Srgrimes if (vq == 0) { 991573Srgrimes /* divide by zero. */ 1001573Srgrimes static volatile const unsigned int zero = 0; 1011573Srgrimes 1021573Srgrimes tmp.ul[H] = tmp.ul[L] = 1 / zero; 1031573Srgrimes if (arq) 1041573Srgrimes *arq = uq; 1051573Srgrimes return (tmp.q); 1061573Srgrimes } 1071573Srgrimes if (uq < vq) { 1081573Srgrimes if (arq) 1091573Srgrimes *arq = uq; 1101573Srgrimes return (0); 1111573Srgrimes } 1121573Srgrimes u = &uspace[0]; 1131573Srgrimes v = &vspace[0]; 1141573Srgrimes q = &qspace[0]; 1151573Srgrimes 1161573Srgrimes /* 1171573Srgrimes * Break dividend and divisor into digits in base B, then 1181573Srgrimes * count leading zeros to determine m and n. When done, we 1191573Srgrimes * will have: 1201573Srgrimes * u = (u[1]u[2]...u[m+n]) sub B 1211573Srgrimes * v = (v[1]v[2]...v[n]) sub B 1221573Srgrimes * v[1] != 0 1231573Srgrimes * 1 < n <= 4 (if n = 1, we use a different division algorithm) 1241573Srgrimes * m >= 0 (otherwise u < v, which we already checked) 1251573Srgrimes * m + n = 4 1261573Srgrimes * and thus 1271573Srgrimes * m = 4 - n <= 2 1281573Srgrimes */ 1291573Srgrimes tmp.uq = uq; 1301573Srgrimes u[0] = 0; 1311573Srgrimes u[1] = HHALF(tmp.ul[H]); 1321573Srgrimes u[2] = LHALF(tmp.ul[H]); 1331573Srgrimes u[3] = HHALF(tmp.ul[L]); 1341573Srgrimes u[4] = LHALF(tmp.ul[L]); 1351573Srgrimes tmp.uq = vq; 1361573Srgrimes v[1] = HHALF(tmp.ul[H]); 1371573Srgrimes v[2] = LHALF(tmp.ul[H]); 1381573Srgrimes v[3] = HHALF(tmp.ul[L]); 1391573Srgrimes v[4] = LHALF(tmp.ul[L]); 1401573Srgrimes for (n = 4; v[1] == 0; v++) { 1411573Srgrimes if (--n == 1) { 1421573Srgrimes u_long rbj; /* r*B+u[j] (not root boy jim) */ 1431573Srgrimes digit q1, q2, q3, q4; 1441573Srgrimes 1451573Srgrimes /* 1461573Srgrimes * Change of plan, per exercise 16. 1471573Srgrimes * r = 0; 1481573Srgrimes * for j = 1..4: 1491573Srgrimes * q[j] = floor((r*B + u[j]) / v), 1501573Srgrimes * r = (r*B + u[j]) % v; 1511573Srgrimes * We unroll this completely here. 1521573Srgrimes */ 1531573Srgrimes t = v[2]; /* nonzero, by definition */ 1541573Srgrimes q1 = u[1] / t; 1551573Srgrimes rbj = COMBINE(u[1] % t, u[2]); 1561573Srgrimes q2 = rbj / t; 1571573Srgrimes rbj = COMBINE(rbj % t, u[3]); 1581573Srgrimes q3 = rbj / t; 1591573Srgrimes rbj = COMBINE(rbj % t, u[4]); 1601573Srgrimes q4 = rbj / t; 1611573Srgrimes if (arq) 1621573Srgrimes *arq = rbj % t; 1631573Srgrimes tmp.ul[H] = COMBINE(q1, q2); 1641573Srgrimes tmp.ul[L] = COMBINE(q3, q4); 1651573Srgrimes return (tmp.q); 1661573Srgrimes } 1671573Srgrimes } 1681573Srgrimes 1691573Srgrimes /* 1701573Srgrimes * By adjusting q once we determine m, we can guarantee that 1711573Srgrimes * there is a complete four-digit quotient at &qspace[1] when 1721573Srgrimes * we finally stop. 1731573Srgrimes */ 1741573Srgrimes for (m = 4 - n; u[1] == 0; u++) 1751573Srgrimes m--; 1761573Srgrimes for (i = 4 - m; --i >= 0;) 1771573Srgrimes q[i] = 0; 1781573Srgrimes q += 4 - m; 1791573Srgrimes 1801573Srgrimes /* 1811573Srgrimes * Here we run Program D, translated from MIX to C and acquiring 1821573Srgrimes * a few minor changes. 1831573Srgrimes * 1841573Srgrimes * D1: choose multiplier 1 << d to ensure v[1] >= B/2. 1851573Srgrimes */ 1861573Srgrimes d = 0; 1871573Srgrimes for (t = v[1]; t < B / 2; t <<= 1) 1881573Srgrimes d++; 1891573Srgrimes if (d > 0) { 1901573Srgrimes shl(&u[0], m + n, d); /* u <<= d */ 1911573Srgrimes shl(&v[1], n - 1, d); /* v <<= d */ 1921573Srgrimes } 1931573Srgrimes /* 1941573Srgrimes * D2: j = 0. 1951573Srgrimes */ 1961573Srgrimes j = 0; 1971573Srgrimes v1 = v[1]; /* for D3 -- note that v[1..n] are constant */ 1981573Srgrimes v2 = v[2]; /* for D3 */ 1991573Srgrimes do { 2001573Srgrimes register digit uj0, uj1, uj2; 2011573Srgrimes 2021573Srgrimes /* 2031573Srgrimes * D3: Calculate qhat (\^q, in TeX notation). 2041573Srgrimes * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and 2051573Srgrimes * let rhat = (u[j]*B + u[j+1]) mod v[1]. 2061573Srgrimes * While rhat < B and v[2]*qhat > rhat*B+u[j+2], 2071573Srgrimes * decrement qhat and increase rhat correspondingly. 2081573Srgrimes * Note that if rhat >= B, v[2]*qhat < rhat*B. 2091573Srgrimes */ 2101573Srgrimes uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */ 2111573Srgrimes uj1 = u[j + 1]; /* for D3 only */ 2121573Srgrimes uj2 = u[j + 2]; /* for D3 only */ 2131573Srgrimes if (uj0 == v1) { 2141573Srgrimes qhat = B; 2151573Srgrimes rhat = uj1; 2161573Srgrimes goto qhat_too_big; 2171573Srgrimes } else { 2181573Srgrimes u_long n = COMBINE(uj0, uj1); 2191573Srgrimes qhat = n / v1; 2201573Srgrimes rhat = n % v1; 2211573Srgrimes } 2221573Srgrimes while (v2 * qhat > COMBINE(rhat, uj2)) { 2231573Srgrimes qhat_too_big: 2241573Srgrimes qhat--; 2251573Srgrimes if ((rhat += v1) >= B) 2261573Srgrimes break; 2271573Srgrimes } 2281573Srgrimes /* 2291573Srgrimes * D4: Multiply and subtract. 2301573Srgrimes * The variable `t' holds any borrows across the loop. 2311573Srgrimes * We split this up so that we do not require v[0] = 0, 2321573Srgrimes * and to eliminate a final special case. 2331573Srgrimes */ 2341573Srgrimes for (t = 0, i = n; i > 0; i--) { 2351573Srgrimes t = u[i + j] - v[i] * qhat - t; 2361573Srgrimes u[i + j] = LHALF(t); 2371573Srgrimes t = (B - HHALF(t)) & (B - 1); 2381573Srgrimes } 2391573Srgrimes t = u[j] - t; 2401573Srgrimes u[j] = LHALF(t); 2411573Srgrimes /* 2421573Srgrimes * D5: test remainder. 2431573Srgrimes * There is a borrow if and only if HHALF(t) is nonzero; 2441573Srgrimes * in that (rare) case, qhat was too large (by exactly 1). 2451573Srgrimes * Fix it by adding v[1..n] to u[j..j+n]. 2461573Srgrimes */ 2471573Srgrimes if (HHALF(t)) { 2481573Srgrimes qhat--; 2491573Srgrimes for (t = 0, i = n; i > 0; i--) { /* D6: add back. */ 2501573Srgrimes t += u[i + j] + v[i]; 2511573Srgrimes u[i + j] = LHALF(t); 2521573Srgrimes t = HHALF(t); 2531573Srgrimes } 2541573Srgrimes u[j] = LHALF(u[j] + t); 2551573Srgrimes } 2561573Srgrimes q[j] = qhat; 2571573Srgrimes } while (++j <= m); /* D7: loop on j. */ 2581573Srgrimes 2591573Srgrimes /* 2601573Srgrimes * If caller wants the remainder, we have to calculate it as 2611573Srgrimes * u[m..m+n] >> d (this is at most n digits and thus fits in 2621573Srgrimes * u[m+1..m+n], but we may need more source digits). 2631573Srgrimes */ 2641573Srgrimes if (arq) { 2651573Srgrimes if (d) { 2661573Srgrimes for (i = m + n; i > m; --i) 2671573Srgrimes u[i] = (u[i] >> d) | 2681573Srgrimes LHALF(u[i - 1] << (HALF_BITS - d)); 2691573Srgrimes u[i] = 0; 2701573Srgrimes } 2711573Srgrimes tmp.ul[H] = COMBINE(uspace[1], uspace[2]); 2721573Srgrimes tmp.ul[L] = COMBINE(uspace[3], uspace[4]); 2731573Srgrimes *arq = tmp.q; 2741573Srgrimes } 2751573Srgrimes 2761573Srgrimes tmp.ul[H] = COMBINE(qspace[1], qspace[2]); 2771573Srgrimes tmp.ul[L] = COMBINE(qspace[3], qspace[4]); 2781573Srgrimes return (tmp.q); 2791573Srgrimes} 280