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