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