Annotation of sys/lib/libkern/qdivrem.c, Revision 1.1
1.1 ! nbrk 1: /*-
! 2: * Copyright (c) 1992, 1993
! 3: * The Regents of the University of California. 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. Neither the name of the University nor the names of its contributors
! 18: * may be used to endorse or promote products derived from this software
! 19: * without specific prior written permission.
! 20: *
! 21: * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
! 22: * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! 23: * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! 24: * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
! 25: * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! 26: * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
! 27: * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
! 28: * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! 29: * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
! 30: * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
! 31: * SUCH DAMAGE.
! 32: */
! 33:
! 34: #if defined(LIBC_SCCS) && !defined(lint)
! 35: static char rcsid[] = "$OpenBSD: qdivrem.c,v 1.8 2005/02/13 03:37:14 jsg Exp $";
! 36: #endif /* LIBC_SCCS and not lint */
! 37:
! 38: /*
! 39: * Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed),
! 40: * section 4.3.1, pp. 257--259.
! 41: */
! 42:
! 43: #include "quad.h"
! 44:
! 45: #define B ((int)1 << HALF_BITS) /* digit base */
! 46:
! 47: /* Combine two `digits' to make a single two-digit number. */
! 48: #define COMBINE(a, b) (((u_int)(a) << HALF_BITS) | (b))
! 49:
! 50: /* select a type for digits in base B: use unsigned short if they fit */
! 51: #if UINT_MAX == 0xffffffffU && USHRT_MAX >= 0xffff
! 52: typedef unsigned short digit;
! 53: #else
! 54: typedef u_int digit;
! 55: #endif
! 56:
! 57: static void shl(digit *p, int len, int sh);
! 58:
! 59: /*
! 60: * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
! 61: *
! 62: * We do this in base 2-sup-HALF_BITS, so that all intermediate products
! 63: * fit within u_int. As a consequence, the maximum length dividend and
! 64: * divisor are 4 `digits' in this base (they are shorter if they have
! 65: * leading zeros).
! 66: */
! 67: u_quad_t
! 68: __qdivrem(u_quad_t uq, u_quad_t vq, u_quad_t *arq)
! 69: {
! 70: union uu tmp;
! 71: digit *u, *v, *q;
! 72: digit v1, v2;
! 73: u_int qhat, rhat, t;
! 74: int m, n, d, j, i;
! 75: digit uspace[5], vspace[5], qspace[5];
! 76:
! 77: /*
! 78: * Take care of special cases: divide by zero, and u < v.
! 79: */
! 80: if (vq == 0) {
! 81: /* divide by zero. */
! 82: static volatile const unsigned int zero = 0;
! 83:
! 84: tmp.ul[H] = tmp.ul[L] = 1 / zero;
! 85: if (arq)
! 86: *arq = uq;
! 87: return (tmp.q);
! 88: }
! 89: if (uq < vq) {
! 90: if (arq)
! 91: *arq = uq;
! 92: return (0);
! 93: }
! 94: u = &uspace[0];
! 95: v = &vspace[0];
! 96: q = &qspace[0];
! 97:
! 98: /*
! 99: * Break dividend and divisor into digits in base B, then
! 100: * count leading zeros to determine m and n. When done, we
! 101: * will have:
! 102: * u = (u[1]u[2]...u[m+n]) sub B
! 103: * v = (v[1]v[2]...v[n]) sub B
! 104: * v[1] != 0
! 105: * 1 < n <= 4 (if n = 1, we use a different division algorithm)
! 106: * m >= 0 (otherwise u < v, which we already checked)
! 107: * m + n = 4
! 108: * and thus
! 109: * m = 4 - n <= 2
! 110: */
! 111: tmp.uq = uq;
! 112: u[0] = 0;
! 113: u[1] = (digit)HHALF(tmp.ul[H]);
! 114: u[2] = (digit)LHALF(tmp.ul[H]);
! 115: u[3] = (digit)HHALF(tmp.ul[L]);
! 116: u[4] = (digit)LHALF(tmp.ul[L]);
! 117: tmp.uq = vq;
! 118: v[1] = (digit)HHALF(tmp.ul[H]);
! 119: v[2] = (digit)LHALF(tmp.ul[H]);
! 120: v[3] = (digit)HHALF(tmp.ul[L]);
! 121: v[4] = (digit)LHALF(tmp.ul[L]);
! 122: for (n = 4; v[1] == 0; v++) {
! 123: if (--n == 1) {
! 124: u_int rbj; /* r*B+u[j] (not root boy jim) */
! 125: digit q1, q2, q3, q4;
! 126:
! 127: /*
! 128: * Change of plan, per exercise 16.
! 129: * r = 0;
! 130: * for j = 1..4:
! 131: * q[j] = floor((r*B + u[j]) / v),
! 132: * r = (r*B + u[j]) % v;
! 133: * We unroll this completely here.
! 134: */
! 135: t = v[2]; /* nonzero, by definition */
! 136: q1 = (digit)(u[1] / t);
! 137: rbj = COMBINE(u[1] % t, u[2]);
! 138: q2 = (digit)(rbj / t);
! 139: rbj = COMBINE(rbj % t, u[3]);
! 140: q3 = (digit)(rbj / t);
! 141: rbj = COMBINE(rbj % t, u[4]);
! 142: q4 = (digit)(rbj / t);
! 143: if (arq)
! 144: *arq = rbj % t;
! 145: tmp.ul[H] = COMBINE(q1, q2);
! 146: tmp.ul[L] = COMBINE(q3, q4);
! 147: return (tmp.q);
! 148: }
! 149: }
! 150:
! 151: /*
! 152: * By adjusting q once we determine m, we can guarantee that
! 153: * there is a complete four-digit quotient at &qspace[1] when
! 154: * we finally stop.
! 155: */
! 156: for (m = 4 - n; u[1] == 0; u++)
! 157: m--;
! 158: for (i = 4 - m; --i >= 0;)
! 159: q[i] = 0;
! 160: q += 4 - m;
! 161:
! 162: /*
! 163: * Here we run Program D, translated from MIX to C and acquiring
! 164: * a few minor changes.
! 165: *
! 166: * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
! 167: */
! 168: d = 0;
! 169: for (t = v[1]; t < B / 2; t <<= 1)
! 170: d++;
! 171: if (d > 0) {
! 172: shl(&u[0], m + n, d); /* u <<= d */
! 173: shl(&v[1], n - 1, d); /* v <<= d */
! 174: }
! 175: /*
! 176: * D2: j = 0.
! 177: */
! 178: j = 0;
! 179: v1 = v[1]; /* for D3 -- note that v[1..n] are constant */
! 180: v2 = v[2]; /* for D3 */
! 181: do {
! 182: digit uj0, uj1, uj2;
! 183:
! 184: /*
! 185: * D3: Calculate qhat (\^q, in TeX notation).
! 186: * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
! 187: * let rhat = (u[j]*B + u[j+1]) mod v[1].
! 188: * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
! 189: * decrement qhat and increase rhat correspondingly.
! 190: * Note that if rhat >= B, v[2]*qhat < rhat*B.
! 191: */
! 192: uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */
! 193: uj1 = u[j + 1]; /* for D3 only */
! 194: uj2 = u[j + 2]; /* for D3 only */
! 195: if (uj0 == v1) {
! 196: qhat = B;
! 197: rhat = uj1;
! 198: goto qhat_too_big;
! 199: } else {
! 200: u_int nn = COMBINE(uj0, uj1);
! 201: qhat = nn / v1;
! 202: rhat = nn % v1;
! 203: }
! 204: while (v2 * qhat > COMBINE(rhat, uj2)) {
! 205: qhat_too_big:
! 206: qhat--;
! 207: if ((rhat += v1) >= B)
! 208: break;
! 209: }
! 210: /*
! 211: * D4: Multiply and subtract.
! 212: * The variable `t' holds any borrows across the loop.
! 213: * We split this up so that we do not require v[0] = 0,
! 214: * and to eliminate a final special case.
! 215: */
! 216: for (t = 0, i = n; i > 0; i--) {
! 217: t = u[i + j] - v[i] * qhat - t;
! 218: u[i + j] = (digit)LHALF(t);
! 219: t = (B - HHALF(t)) & (B - 1);
! 220: }
! 221: t = u[j] - t;
! 222: u[j] = (digit)LHALF(t);
! 223: /*
! 224: * D5: test remainder.
! 225: * There is a borrow if and only if HHALF(t) is nonzero;
! 226: * in that (rare) case, qhat was too large (by exactly 1).
! 227: * Fix it by adding v[1..n] to u[j..j+n].
! 228: */
! 229: if (HHALF(t)) {
! 230: qhat--;
! 231: for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
! 232: t += u[i + j] + v[i];
! 233: u[i + j] = (digit)LHALF(t);
! 234: t = HHALF(t);
! 235: }
! 236: u[j] = (digit)LHALF(u[j] + t);
! 237: }
! 238: q[j] = (digit)qhat;
! 239: } while (++j <= m); /* D7: loop on j. */
! 240:
! 241: /*
! 242: * If caller wants the remainder, we have to calculate it as
! 243: * u[m..m+n] >> d (this is at most n digits and thus fits in
! 244: * u[m+1..m+n], but we may need more source digits).
! 245: */
! 246: if (arq) {
! 247: if (d) {
! 248: for (i = m + n; i > m; --i)
! 249: u[i] = (digit)(((u_int)u[i] >> d) |
! 250: LHALF((u_int)u[i - 1] << (HALF_BITS - d)));
! 251: u[i] = 0;
! 252: }
! 253: tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
! 254: tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
! 255: *arq = tmp.q;
! 256: }
! 257:
! 258: tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
! 259: tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
! 260: return (tmp.q);
! 261: }
! 262:
! 263: /*
! 264: * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
! 265: * `fall out' the left (there never will be any such anyway).
! 266: * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS.
! 267: */
! 268: static void
! 269: shl(digit *p, int len, int sh)
! 270: {
! 271: int i;
! 272:
! 273: for (i = 0; i < len; i++)
! 274: p[i] = (digit)(LHALF((u_int)p[i] << sh) |
! 275: ((u_int)p[i + 1] >> (HALF_BITS - sh)));
! 276: p[i] = (digit)(LHALF((u_int)p[i] << sh));
! 277: }
CVSweb
![[BACK]](/icons/back.gif){kind=icon}
Return to qdivrem.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / lib / libkern