qdivrem.c

Version: - FREEBSD - FREEBSD-13-STABLE - FREEBSD-13-0 - FREEBSD-12-STABLE - FREEBSD-12-0 - FREEBSD-11-STABLE - FREEBSD-11-0 - FREEBSD-10-STABLE - FREEBSD-10-0 - FREEBSD-9-STABLE - FREEBSD-9-0 - FREEBSD-8-STABLE - FREEBSD-8-0 - FREEBSD-7-STABLE - FREEBSD-7-0 - FREEBSD-6-STABLE - FREEBSD-6-0 - FREEBSD-5-STABLE - FREEBSD-5-0 - FREEBSD-4-STABLE - FREEBSD-3-STABLE - FREEBSD22 - l41 - OPENBSD - linux-2.6 - MK84 - PLAN9 - xnu-8792
SearchContext: - none - 3 - 10
    1 /*      $NetBSD: qdivrem.c,v 1.10 2003/08/07 16:32:10 agc Exp $ */
    2 
    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 #include <sys/cdefs.h>
   37 #if defined(LIBC_SCCS) && !defined(lint)
   38 #if 0
   39 static char sccsid[] = "@(#)qdivrem.c   8.1 (Berkeley) 6/4/93";
   40 #else
   41 __RCSID("$NetBSD: qdivrem.c,v 1.10 2003/08/07 16:32:10 agc Exp $");
   42 #endif
   43 #endif /* LIBC_SCCS and not lint */
   44 
   45 /*
   46  * Multiprecision divide.  This algorithm is from Knuth vol. 2 (2nd ed),
   47  * section 4.3.1, pp. 257--259.
   48  */
   49 
   50 #include "quad.h"
   51 
   52 #define B       ((int)1 << HALF_BITS)   /* digit base */
   53 
   54 /* Combine two `digits' to make a single two-digit number. */
   55 #define COMBINE(a, b) (((u_int)(a) << HALF_BITS) | (b))
   56 
   57 /* select a type for digits in base B: use unsigned short if they fit */
   58 #if UINT_MAX == 0xffffffffU && USHRT_MAX >= 0xffff
   59 typedef unsigned short digit;
   60 #else
   61 typedef u_int digit;
   62 #endif
   63 
   64 static void shl __P((digit *p, int len, int sh));
   65 
   66 /*
   67  * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
   68  *
   69  * We do this in base 2-sup-HALF_BITS, so that all intermediate products
   70  * fit within u_int.  As a consequence, the maximum length dividend and
   71  * divisor are 4 `digits' in this base (they are shorter if they have
   72  * leading zeros).
   73  */
   74 u_quad_t
   75 __qdivrem(uq, vq, arq)
   76         u_quad_t uq, vq, *arq;
   77 {
   78         union uu tmp;
   79         digit *u, *v, *q;
   80         digit v1, v2;
   81         u_int qhat, rhat, t;
   82         int m, n, d, j, i;
   83         digit uspace[5], vspace[5], qspace[5];
   84 
   85         /*
   86          * Take care of special cases: divide by zero, and u < v.
   87          */
   88         if (vq == 0) {
   89                 /* divide by zero. */
   90                 static volatile const unsigned int zero = 0;
   91 
   92                 tmp.ul[H] = tmp.ul[L] = 1 / zero;
   93                 if (arq)
   94                         *arq = uq;
   95                 return (tmp.q);
   96         }
   97         if (uq < vq) {
   98                 if (arq)
   99                         *arq = uq;
  100                 return (0);
  101         }
  102         u = &uspace[0];
  103         v = &vspace[0];
  104         q = &qspace[0];
  105 
  106         /*
  107          * Break dividend and divisor into digits in base B, then
  108          * count leading zeros to determine m and n.  When done, we
  109          * will have:
  110          *      u = (u[1]u[2]...u[m+n]) sub B
  111          *      v = (v[1]v[2]...v[n]) sub B
  112          *      v[1] != 0
  113          *      1 < n <= 4 (if n = 1, we use a different division algorithm)
  114          *      m >= 0 (otherwise u < v, which we already checked)
  115          *      m + n = 4
  116          * and thus
  117          *      m = 4 - n <= 2
  118          */
  119         tmp.uq = uq;
  120         u[0] = 0;
  121         u[1] = (digit)HHALF(tmp.ul[H]);
  122         u[2] = (digit)LHALF(tmp.ul[H]);
  123         u[3] = (digit)HHALF(tmp.ul[L]);
  124         u[4] = (digit)LHALF(tmp.ul[L]);
  125         tmp.uq = vq;
  126         v[1] = (digit)HHALF(tmp.ul[H]);
  127         v[2] = (digit)LHALF(tmp.ul[H]);
  128         v[3] = (digit)HHALF(tmp.ul[L]);
  129         v[4] = (digit)LHALF(tmp.ul[L]);
  130         for (n = 4; v[1] == 0; v++) {
  131                 if (--n == 1) {
  132                         u_int rbj;      /* r*B+u[j] (not root boy jim) */
  133                         digit q1, q2, q3, q4;
  134 
  135                         /*
  136                          * Change of plan, per exercise 16.
  137                          *      r = 0;
  138                          *      for j = 1..4:
  139                          *              q[j] = floor((r*B + u[j]) / v),
  140                          *              r = (r*B + u[j]) % v;
  141                          * We unroll this completely here.
  142                          */
  143                         t = v[2];       /* nonzero, by definition */
  144                         q1 = (digit)(u[1] / t);
  145                         rbj = COMBINE(u[1] % t, u[2]);
  146                         q2 = (digit)(rbj / t);
  147                         rbj = COMBINE(rbj % t, u[3]);
  148                         q3 = (digit)(rbj / t);
  149                         rbj = COMBINE(rbj % t, u[4]);
  150                         q4 = (digit)(rbj / t);
  151                         if (arq)
  152                                 *arq = rbj % t;
  153                         tmp.ul[H] = COMBINE(q1, q2);
  154                         tmp.ul[L] = COMBINE(q3, q4);
  155                         return (tmp.q);
  156                 }
  157         }
  158 
  159         /*
  160          * By adjusting q once we determine m, we can guarantee that
  161          * there is a complete four-digit quotient at &qspace[1] when
  162          * we finally stop.
  163          */
  164         for (m = 4 - n; u[1] == 0; u++)
  165                 m--;
  166         for (i = 4 - m; --i >= 0;)
  167                 q[i] = 0;
  168         q += 4 - m;
  169 
  170         /*
  171          * Here we run Program D, translated from MIX to C and acquiring
  172          * a few minor changes.
  173          *
  174          * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
  175          */
  176         d = 0;
  177         for (t = v[1]; t < B / 2; t <<= 1)
  178                 d++;
  179         if (d > 0) {
  180                 shl(&u[0], m + n, d);           /* u <<= d */
  181                 shl(&v[1], n - 1, d);           /* v <<= d */
  182         }
  183         /*
  184          * D2: j = 0.
  185          */
  186         j = 0;
  187         v1 = v[1];      /* for D3 -- note that v[1..n] are constant */
  188         v2 = v[2];      /* for D3 */
  189         do {
  190                 digit uj0, uj1, uj2;
  191                 
  192                 /*
  193                  * D3: Calculate qhat (\^q, in TeX notation).
  194                  * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
  195                  * let rhat = (u[j]*B + u[j+1]) mod v[1].
  196                  * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
  197                  * decrement qhat and increase rhat correspondingly.
  198                  * Note that if rhat >= B, v[2]*qhat < rhat*B.
  199                  */
  200                 uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */
  201                 uj1 = u[j + 1]; /* for D3 only */
  202                 uj2 = u[j + 2]; /* for D3 only */
  203                 if (uj0 == v1) {
  204                         qhat = B;
  205                         rhat = uj1;
  206                         goto qhat_too_big;
  207                 } else {
  208                         u_int nn = COMBINE(uj0, uj1);
  209                         qhat = nn / v1;
  210                         rhat = nn % v1;
  211                 }
  212                 while (v2 * qhat > COMBINE(rhat, uj2)) {
  213         qhat_too_big:
  214                         qhat--;
  215                         if ((rhat += v1) >= B)
  216                                 break;
  217                 }
  218                 /*
  219                  * D4: Multiply and subtract.
  220                  * The variable `t' holds any borrows across the loop.
  221                  * We split this up so that we do not require v[0] = 0,
  222                  * and to eliminate a final special case.
  223                  */
  224                 for (t = 0, i = n; i > 0; i--) {
  225                         t = u[i + j] - v[i] * qhat - t;
  226                         u[i + j] = (digit)LHALF(t);
  227                         t = (B - HHALF(t)) & (B - 1);
  228                 }
  229                 t = u[j] - t;
  230                 u[j] = (digit)LHALF(t);
  231                 /*
  232                  * D5: test remainder.
  233                  * There is a borrow if and only if HHALF(t) is nonzero;
  234                  * in that (rare) case, qhat was too large (by exactly 1).
  235                  * Fix it by adding v[1..n] to u[j..j+n].
  236                  */
  237                 if (HHALF(t)) {
  238                         qhat--;
  239                         for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
  240                                 t += u[i + j] + v[i];
  241                                 u[i + j] = (digit)LHALF(t);
  242                                 t = HHALF(t);
  243                         }
  244                         u[j] = (digit)LHALF(u[j] + t);
  245                 }
  246                 q[j] = (digit)qhat;
  247         } while (++j <= m);             /* D7: loop on j. */
  248 
  249         /*
  250          * If caller wants the remainder, we have to calculate it as
  251          * u[m..m+n] >> d (this is at most n digits and thus fits in
  252          * u[m+1..m+n], but we may need more source digits).
  253          */
  254         if (arq) {
  255                 if (d) {
  256                         for (i = m + n; i > m; --i)
  257                                 u[i] = (digit)(((u_int)u[i] >> d) |
  258                                     LHALF((u_int)u[i - 1] << (HALF_BITS - d)));
  259                         u[i] = 0;
  260                 }
  261                 tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
  262                 tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
  263                 *arq = tmp.q;
  264         }
  265 
  266         tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
  267         tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
  268         return (tmp.q);
  269 }
  270 
  271 /*
  272  * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
  273  * `fall out' the left (there never will be any such anyway).
  274  * We may assume len >= 0.  NOTE THAT THIS WRITES len+1 DIGITS.
  275  */
  276 static void
  277 shl(digit *p, int len, int sh)
  278 {
  279         int i;
  280 
  281         for (i = 0; i < len; i++)
  282                 p[i] = (digit)(LHALF((u_int)p[i] << sh) |
  283                     ((u_int)p[i + 1] >> (HALF_BITS - sh)));
  284         p[i] = (digit)(LHALF((u_int)p[i] << sh));
  285 }
Cache object: b9b37af5522de37bb5ded8b0c861e484
[ source navigation ] [ diff markup ] [ identifier search ] [ freetext search ] [ file search ] [ list types ] [ track identifier ]
This page is part of the FreeBSD/Linux Linux Kernel Cross-Reference, and was automatically generated using a modified version of the LXR engine.
FreeBSD/Linux Kernel Cross Reference sys/lib/libkern/qdivrem.c

FreeBSD/Linux Kernel Cross Reference
sys/lib/libkern/qdivrem.c
