Subversion Repositories shark

/*-

 * Copyright (c) 1992, 1993

 *      The Regents of the University of California.  All rights reserved.

 *

 * This software was developed by the Computer Systems Engineering group

 * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and

 * contributed to Berkeley.

 *

 * Redistribution and use in source and binary forms, with or without

 * modification, are permitted provided that the following conditions

 * are met:

 * 1. Redistributions of source code must retain the above copyright

 *    notice, this list of conditions and the following disclaimer.

 * 2. Redistributions in binary form must reproduce the above copyright

 *    notice, this list of conditions and the following disclaimer in the

 *    documentation and/or other materials provided with the distribution.

 * 3. All advertising materials mentioning features or use of this software

 *    must display the following acknowledgement:

 *      This product includes software developed by the University of

 *      California, Berkeley and its contributors.

 * 4. Neither the name of the University nor the names of its contributors

 *    may be used to endorse or promote products derived from this software

 *    without specific prior written permission.

 *

 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND

 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE

 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS

 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY

 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF

 * SUCH DAMAGE.

 */

#if defined(LIBC_SCCS) && !defined(lint)

static char sccsid[] = "@(#)qdivrem.c   8.1 (Berkeley) 6/4/93";

#endif /* LIBC_SCCS and not lint */

/*

 * Multiprecision divide.  This algorithm is from Knuth vol. 2 (2nd ed),

 * section 4.3.1, pp. 257--259.

 */

#include "quad.h"

#define B       (1 << HALF_BITS)        /* digit base */

/* Combine two `digits' to make a single two-digit number. */

#define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b))

/* select a type for digits in base B: use unsigned short if they fit */

#if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff

typedef unsigned short digit;

#else

typedef u_long digit;

#endif

/*

 * Shift p[0]..p[len] left `sh' bits, ignoring any bits that

 * `fall out' the left (there never will be any such anyway).

 * We may assume len >= 0.  NOTE THAT THIS WRITES len+1 DIGITS.

 */

static void

shl(register digit *p, register int len, register int sh)

{

        register int i;

        for (i = 0; i < len; i++)

                p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh));

        p[i] = LHALF(p[i] << sh);

}

/*

 * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.

 *

 * We do this in base 2-sup-HALF_BITS, so that all intermediate products

 * fit within u_long.  As a consequence, the maximum length dividend and

 * divisor are 4 `digits' in this base (they are shorter if they have

 * leading zeros).

 */

u_quad_t

__qdivrem(uq, vq, arq)

        u_quad_t uq, vq, *arq;

{

        union uu tmp;

        digit *u, *v, *q;

        register digit v1, v2;

        u_long qhat, rhat, t;

        int m, n, d, j, i;

        digit uspace[5], vspace[5], qspace[5];

        /*

         * Take care of special cases: divide by zero, and u < v.

         */

        if (vq == 0) {

                /* divide by zero. */

                static volatile const unsigned int zero = 0;

                tmp.ul[H] = tmp.ul[L] = 1 / zero;

                if (arq)

                        *arq = uq;

                return (tmp.q);

        }

        if (uq < vq) {

                if (arq)

                        *arq = uq;

                return (0);

        }

        u = &uspace[0];

        v = &vspace[0];

        q = &qspace[0];

        /*

         * Break dividend and divisor into digits in base B, then

         * count leading zeros to determine m and n.  When done, we

         * will have:

         *      u = (u[1]u[2]...u[m+n]) sub B

         *      v = (v[1]v[2]...v[n]) sub B

         *      v[1] != 0

         *      1 < n <= 4 (if n = 1, we use a different division algorithm)

         *      m >= 0 (otherwise u < v, which we already checked)

         *      m + n = 4

         * and thus

         *      m = 4 - n <= 2

         */

        tmp.uq = uq;

        u[0] = 0;

        u[1] = HHALF(tmp.ul[H]);

        u[2] = LHALF(tmp.ul[H]);

        u[3] = HHALF(tmp.ul[L]);

        u[4] = LHALF(tmp.ul[L]);

        tmp.uq = vq;

        v[1] = HHALF(tmp.ul[H]);

        v[2] = LHALF(tmp.ul[H]);

        v[3] = HHALF(tmp.ul[L]);

        v[4] = LHALF(tmp.ul[L]);

        for (n = 4; v[1] == 0; v++) {

                if (--n == 1) {

                        u_long rbj;     /* r*B+u[j] (not root boy jim) */

                        digit q1, q2, q3, q4;

                        /*

                         * Change of plan, per exercise 16.

                         *      r = 0;

                         *      for j = 1..4:

                         *              q[j] = floor((r*B + u[j]) / v),

                         *              r = (r*B + u[j]) % v;

                         * We unroll this completely here.

                         */

                        t = v[2];       /* nonzero, by definition */

                        q1 = u[1] / t;

                        rbj = COMBINE(u[1] % t, u[2]);

                        q2 = rbj / t;

                        rbj = COMBINE(rbj % t, u[3]);

                        q3 = rbj / t;

                        rbj = COMBINE(rbj % t, u[4]);

                        q4 = rbj / t;

                        if (arq)

                                *arq = rbj % t;

                        tmp.ul[H] = COMBINE(q1, q2);

                        tmp.ul[L] = COMBINE(q3, q4);

                        return (tmp.q);

                }

        }

        /*

         * By adjusting q once we determine m, we can guarantee that

         * there is a complete four-digit quotient at &qspace[1] when

         * we finally stop.

         */

        for (m = 4 - n; u[1] == 0; u++)

                m--;

        for (i = 4 - m; --i >= 0;)

                q[i] = 0;

        q += 4 - m;

        /*

         * Here we run Program D, translated from MIX to C and acquiring

         * a few minor changes.

         *

         * D1: choose multiplier 1 << d to ensure v[1] >= B/2.

         */

        d = 0;

        for (t = v[1]; t < B / 2; t <<= 1)

                d++;

        if (d > 0) {

                shl(&u[0], m + n, d);           /* u <<= d */

                shl(&v[1], n - 1, d);           /* v <<= d */

        }

        /*

         * D2: j = 0.

         */

        j = 0;

        v1 = v[1];      /* for D3 -- note that v[1..n] are constant */

        v2 = v[2];      /* for D3 */

        do {

                register digit uj0, uj1, uj2;

                /*

                 * D3: Calculate qhat (\^q, in TeX notation).

                 * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and

                 * let rhat = (u[j]*B + u[j+1]) mod v[1].

                 * While rhat < B and v[2]*qhat > rhat*B+u[j+2],

                 * decrement qhat and increase rhat correspondingly.

                 * Note that if rhat >= B, v[2]*qhat < rhat*B.

                 */

                uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */

                uj1 = u[j + 1]; /* for D3 only */

                uj2 = u[j + 2]; /* for D3 only */

                if (uj0 == v1) {

                        qhat = B;

                        rhat = uj1;

                        goto qhat_too_big;

                } else {

                        u_long n = COMBINE(uj0, uj1);

                        qhat = n / v1;

                        rhat = n % v1;

                }

                while (v2 * qhat > COMBINE(rhat, uj2)) {

        qhat_too_big:

                        qhat--;

                        if ((rhat += v1) >= B)

                                break;

                }

                /*

                 * D4: Multiply and subtract.

                 * The variable `t' holds any borrows across the loop.

                 * We split this up so that we do not require v[0] = 0,

                 * and to eliminate a final special case.

                 */

                for (t = 0, i = n; i > 0; i--) {

                        t = u[i + j] - v[i] * qhat - t;

                        u[i + j] = LHALF(t);

                        t = (B - HHALF(t)) & (B - 1);

                }

                t = u[j] - t;

                u[j] = LHALF(t);

                /*

                 * D5: test remainder.

                 * There is a borrow if and only if HHALF(t) is nonzero;

                 * in that (rare) case, qhat was too large (by exactly 1).

                 * Fix it by adding v[1..n] to u[j..j+n].

                 */

                if (HHALF(t)) {

                        qhat--;

                        for (t = 0, i = n; i > 0; i--) { /* D6: add back. */

                                t += u[i + j] + v[i];

                                u[i + j] = LHALF(t);

                                t = HHALF(t);

                        }

                        u[j] = LHALF(u[j] + t);

                }

                q[j] = qhat;

        } while (++j <= m);             /* D7: loop on j. */

        /*

         * If caller wants the remainder, we have to calculate it as

         * u[m..m+n] >> d (this is at most n digits and thus fits in

         * u[m+1..m+n], but we may need more source digits).

         */

        if (arq) {

                if (d) {

                        for (i = m + n; i > m; --i)

                                u[i] = (u[i] >> d) |

                                    LHALF(u[i - 1] << (HALF_BITS - d));

                        u[i] = 0;

                }

                tmp.ul[H] = COMBINE(uspace[1], uspace[2]);

                tmp.ul[L] = COMBINE(uspace[3], uspace[4]);

                *arq = tmp.q;

        }

        tmp.ul[H] = COMBINE(qspace[1], qspace[2]);

        tmp.ul[L] = COMBINE(qspace[3], qspace[4]);

        return (tmp.q);

}