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.4 2003/06/02 20:18:36 millert 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 ((long)1 << HALF_BITS) /* digit base */ 46 47 /* Combine two `digits' to make a single two-digit number. */ 48 #define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b)) 49 50 /* select a type for digits in base B: use unsigned short if they fit */ 51 #if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff 52 typedef unsigned short digit; 53 #else 54 typedef u_long 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_long. 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(uq, vq, arq) 69 u_quad_t uq, vq, *arq; 70 { 71 union uu tmp; 72 digit *u, *v, *q; 73 register digit v1, v2; 74 u_long qhat, rhat, t; 75 int m, n, d, j, i; 76 digit uspace[5], vspace[5], qspace[5]; 77 78 /* 79 * Take care of special cases: divide by zero, and u < v. 80 */ 81 if (vq == 0) { 82 /* divide by zero. */ 83 static volatile const unsigned int zero = 0; 84 85 tmp.ul[H] = tmp.ul[L] = 1 / zero; 86 if (arq) 87 *arq = uq; 88 return (tmp.q); 89 } 90 if (uq < vq) { 91 if (arq) 92 *arq = uq; 93 return (0); 94 } 95 u = &uspace[0]; 96 v = &vspace[0]; 97 q = &qspace[0]; 98 99 /* 100 * Break dividend and divisor into digits in base B, then 101 * count leading zeros to determine m and n. When done, we 102 * will have: 103 * u = (u[1]u[2]...u[m+n]) sub B 104 * v = (v[1]v[2]...v[n]) sub B 105 * v[1] != 0 106 * 1 < n <= 4 (if n = 1, we use a different division algorithm) 107 * m >= 0 (otherwise u < v, which we already checked) 108 * m + n = 4 109 * and thus 110 * m = 4 - n <= 2 111 */ 112 tmp.uq = uq; 113 u[0] = 0; 114 u[1] = HHALF(tmp.ul[H]); 115 u[2] = LHALF(tmp.ul[H]); 116 u[3] = HHALF(tmp.ul[L]); 117 u[4] = LHALF(tmp.ul[L]); 118 tmp.uq = vq; 119 v[1] = HHALF(tmp.ul[H]); 120 v[2] = LHALF(tmp.ul[H]); 121 v[3] = HHALF(tmp.ul[L]); 122 v[4] = LHALF(tmp.ul[L]); 123 for (n = 4; v[1] == 0; v++) { 124 if (--n == 1) { 125 u_long rbj; /* r*B+u[j] (not root boy jim) */ 126 digit q1, q2, q3, q4; 127 128 /* 129 * Change of plan, per exercise 16. 130 * r = 0; 131 * for j = 1..4: 132 * q[j] = floor((r*B + u[j]) / v), 133 * r = (r*B + u[j]) % v; 134 * We unroll this completely here. 135 */ 136 t = v[2]; /* nonzero, by definition */ 137 q1 = u[1] / t; 138 rbj = COMBINE(u[1] % t, u[2]); 139 q2 = rbj / t; 140 rbj = COMBINE(rbj % t, u[3]); 141 q3 = rbj / t; 142 rbj = COMBINE(rbj % t, u[4]); 143 q4 = rbj / t; 144 if (arq) 145 *arq = rbj % t; 146 tmp.ul[H] = COMBINE(q1, q2); 147 tmp.ul[L] = COMBINE(q3, q4); 148 return (tmp.q); 149 } 150 } 151 152 /* 153 * By adjusting q once we determine m, we can guarantee that 154 * there is a complete four-digit quotient at &qspace[1] when 155 * we finally stop. 156 */ 157 for (m = 4 - n; u[1] == 0; u++) 158 m--; 159 for (i = 4 - m; --i >= 0;) 160 q[i] = 0; 161 q += 4 - m; 162 163 /* 164 * Here we run Program D, translated from MIX to C and acquiring 165 * a few minor changes. 166 * 167 * D1: choose multiplier 1 << d to ensure v[1] >= B/2. 168 */ 169 d = 0; 170 for (t = v[1]; t < B / 2; t <<= 1) 171 d++; 172 if (d > 0) { 173 shl(&u[0], m + n, d); /* u <<= d */ 174 shl(&v[1], n - 1, d); /* v <<= d */ 175 } 176 /* 177 * D2: j = 0. 178 */ 179 j = 0; 180 v1 = v[1]; /* for D3 -- note that v[1..n] are constant */ 181 v2 = v[2]; /* for D3 */ 182 do { 183 register digit uj0, uj1, uj2; 184 185 /* 186 * D3: Calculate qhat (\^q, in TeX notation). 187 * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and 188 * let rhat = (u[j]*B + u[j+1]) mod v[1]. 189 * While rhat < B and v[2]*qhat > rhat*B+u[j+2], 190 * decrement qhat and increase rhat correspondingly. 191 * Note that if rhat >= B, v[2]*qhat < rhat*B. 192 */ 193 uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */ 194 uj1 = u[j + 1]; /* for D3 only */ 195 uj2 = u[j + 2]; /* for D3 only */ 196 if (uj0 == v1) { 197 qhat = B; 198 rhat = uj1; 199 goto qhat_too_big; 200 } else { 201 u_long n = COMBINE(uj0, uj1); 202 qhat = n / v1; 203 rhat = n % v1; 204 } 205 while (v2 * qhat > COMBINE(rhat, uj2)) { 206 qhat_too_big: 207 qhat--; 208 if ((rhat += v1) >= B) 209 break; 210 } 211 /* 212 * D4: Multiply and subtract. 213 * The variable `t' holds any borrows across the loop. 214 * We split this up so that we do not require v[0] = 0, 215 * and to eliminate a final special case. 216 */ 217 for (t = 0, i = n; i > 0; i--) { 218 t = u[i + j] - v[i] * qhat - t; 219 u[i + j] = LHALF(t); 220 t = (B - HHALF(t)) & (B - 1); 221 } 222 t = u[j] - t; 223 u[j] = LHALF(t); 224 /* 225 * D5: test remainder. 226 * There is a borrow if and only if HHALF(t) is nonzero; 227 * in that (rare) case, qhat was too large (by exactly 1). 228 * Fix it by adding v[1..n] to u[j..j+n]. 229 */ 230 if (HHALF(t)) { 231 qhat--; 232 for (t = 0, i = n; i > 0; i--) { /* D6: add back. */ 233 t += u[i + j] + v[i]; 234 u[i + j] = LHALF(t); 235 t = HHALF(t); 236 } 237 u[j] = LHALF(u[j] + t); 238 } 239 q[j] = qhat; 240 } while (++j <= m); /* D7: loop on j. */ 241 242 /* 243 * If caller wants the remainder, we have to calculate it as 244 * u[m..m+n] >> d (this is at most n digits and thus fits in 245 * u[m+1..m+n], but we may need more source digits). 246 */ 247 if (arq) { 248 if (d) { 249 for (i = m + n; i > m; --i) 250 u[i] = (u[i] >> d) | 251 LHALF(u[i - 1] << (HALF_BITS - d)); 252 u[i] = 0; 253 } 254 tmp.ul[H] = COMBINE(uspace[1], uspace[2]); 255 tmp.ul[L] = COMBINE(uspace[3], uspace[4]); 256 *arq = tmp.q; 257 } 258 259 tmp.ul[H] = COMBINE(qspace[1], qspace[2]); 260 tmp.ul[L] = COMBINE(qspace[3], qspace[4]); 261 return (tmp.q); 262 } 263 264 /* 265 * Shift p[0]..p[len] left `sh' bits, ignoring any bits that 266 * `fall out' the left (there never will be any such anyway). 267 * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS. 268 */ 269 static void 270 shl(register digit *p, register int len, register int sh) 271 { 272 register int i; 273 274 for (i = 0; i < len; i++) 275 p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh)); 276 p[i] = LHALF(p[i] << sh); 277 } 278