1 /*- 2 * Copyright (c) 1992 The Regents of the University of California. 3 * All rights reserved. 4 * 5 * %sccs.include.redist.c% 6 */ 7 8 #if defined(LIBC_SCCS) && !defined(lint) 9 static char sccsid[] = "@(#)qdivrem.c 5.3 (Berkeley) 05/12/92"; 10 #endif /* LIBC_SCCS and not lint */ 11 12 #include "longlong.h" 13 14 static int bshift (); 15 16 /* Divide a by b, producing quotient q and remainder r. 17 18 sizeof a is m 19 sizeof b is n 20 sizeof q is m - n 21 sizeof r is n 22 23 The quotient must fit in m - n bytes, i.e., the most significant 24 n digits of a must be less than b, and m must be greater than n. */ 25 26 /* The name of this used to be __div_internal, 27 but that is too long for SYSV. */ 28 29 void 30 __bdiv (a, b, q, r, m, n) 31 unsigned short *a, *b, *q, *r; 32 size_t m, n; 33 { 34 unsigned long qhat, rhat; 35 unsigned long acc; 36 unsigned short *u = (unsigned short *) alloca (m); 37 unsigned short *v = (unsigned short *) alloca (n); 38 unsigned short *u0, *u1, *u2; 39 unsigned short *v0; 40 int d, qn; 41 int i, j; 42 43 m /= sizeof *a; 44 n /= sizeof *b; 45 qn = m - n; 46 47 /* Remove leading zero digits from divisor, and the same number of 48 digits (which must be zero) from dividend. */ 49 50 while (b[big_end (n)] == 0) 51 { 52 r[big_end (n)] = 0; 53 54 a += little_end (2); 55 b += little_end (2); 56 r += little_end (2); 57 m--; 58 n--; 59 60 /* Check for zero divisor. */ 61 if (n == 0) 62 { 63 *r /= n; /* force a divide-by-zero trap */ 64 return; 65 } 66 } 67 68 /* If divisor is a single digit, do short division. */ 69 70 if (n == 1) 71 { 72 acc = a[big_end (m)]; 73 a += little_end (2); 74 for (j = big_end (qn); is_not_lsd (j, qn); j = next_lsd (j)) 75 { 76 acc = (acc << 16) | a[j]; 77 q[j] = acc / *b; 78 acc = acc % *b; 79 } 80 *r = acc; 81 return; 82 } 83 84 /* No such luck, must do long division. Shift divisor and dividend 85 left until the high bit of the divisor is 1. */ 86 87 for (d = 0; d < 16; d++) 88 if (b[big_end (n)] & (1 << (16 - 1 - d))) 89 break; 90 91 bshift (a, d, u, 0, m); 92 bshift (b, d, v, 0, n); 93 94 /* Get pointers to the high dividend and divisor digits. */ 95 96 u0 = u + big_end (m) - big_end (qn); 97 u1 = next_lsd (u0); 98 u2 = next_lsd (u1); 99 u += little_end (2); 100 101 v0 = v + big_end (n); 102 103 /* Main loop: find a quotient digit, multiply it by the divisor, 104 and subtract that from the dividend, shifted over the right amount. */ 105 106 for (j = big_end (qn); is_not_lsd (j, qn); j = next_lsd (j)) 107 { 108 /* Quotient digit initial guess: high 2 dividend digits over high 109 divisor digit. */ 110 111 if (u0[j] == *v0) 112 { 113 qhat = B - 1; 114 rhat = (unsigned long) *v0 + u1[j]; 115 } 116 else 117 { 118 unsigned long numerator = ((unsigned long) u0[j] << 16) | u1[j]; 119 qhat = numerator / *v0; 120 rhat = numerator % *v0; 121 } 122 123 /* Now get the quotient right for high 3 dividend digits over 124 high 2 divisor digits. */ 125 126 while (rhat < B && qhat * *next_lsd (v0) > ((rhat << 16) | u2[j])) 127 { 128 qhat -= 1; 129 rhat += *v0; 130 } 131 132 /* Multiply quotient by divisor, subtract from dividend. */ 133 134 acc = 0; 135 for (i = little_end (n); is_not_msd (i, n); i = next_msd (i)) 136 { 137 acc += (unsigned long) (u + j)[i] - v[i] * qhat; 138 (u + j)[i] = acc & low16; 139 if (acc < B) 140 acc = 0; 141 else 142 acc = (acc >> 16) | -B; 143 } 144 145 q[j] = qhat; 146 147 /* Quotient may have been too high by 1. If dividend went negative, 148 decrement the quotient by 1 and add the divisor back. */ 149 150 if ((signed long) (acc + u0[j]) < 0) 151 { 152 q[j] -= 1; 153 acc = 0; 154 for (i = little_end (n); is_not_msd (i, n); i = next_msd (i)) 155 { 156 acc += (unsigned long) (u + j)[i] + v[i]; 157 (u + j)[i] = acc & low16; 158 acc = acc >> 16; 159 } 160 } 161 } 162 163 /* Now the remainder is what's left of the dividend, shifted right 164 by the amount of the normalizing left shift at the top. */ 165 166 r[big_end (n)] = bshift (u + 1 + little_end (j - 1), 167 16 - d, 168 r + little_end (2), 169 u[little_end (m - 1)] >> d, 170 n - 1); 171 } 172 173 /* Left shift U by K giving W; fill the introduced low-order bits with 174 CARRY_IN. Length of U and W is N. Return carry out. K must be 175 in 0 .. 16. */ 176 177 static int 178 bshift (u, k, w, carry_in, n) 179 unsigned short *u, *w, carry_in; 180 int k, n; 181 { 182 unsigned long acc; 183 int i; 184 185 if (k == 0) 186 { 187 while (n-- > 0) 188 *w++ = *u++; 189 return 0; 190 } 191 192 acc = carry_in; 193 for (i = little_end (n); is_not_msd (i, n); i = next_msd (i)) 194 { 195 acc |= (unsigned long) u[i] << k; 196 w[i] = acc & low16; 197 acc = acc >> 16; 198 } 199 return acc; 200 } 201