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