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.5 (Berkeley) 05/12/92"; 10 #endif /* LIBC_SCCS and not lint */ 11 12 /* Copyright (C) 1989, 1992 Free Software Foundation, Inc. 13 14 This file is part of GNU CC. 15 16 GNU CC is free software; you can redistribute it and/or modify 17 it under the terms of the GNU General Public License as published by 18 the Free Software Foundation; either version 2, or (at your option) 19 any later version. 20 21 GNU CC is distributed in the hope that it will be useful, 22 but WITHOUT ANY WARRANTY; without even the implied warranty of 23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 24 GNU General Public License for more details. 25 26 You should have received a copy of the GNU General Public License 27 along with GNU CC; see the file COPYING. If not, write to 28 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */ 29 30 /* As a special exception, if you link this library with files 31 compiled with GCC to produce an executable, this does not cause 32 the resulting executable to be covered by the GNU General Public License. 33 This exception does not however invalidate any other reasons why 34 the executable file might be covered by the GNU General Public License. */ 35 36 #include <unistd.h> 37 38 #include "longlong.h" 39 40 static int bshift (); 41 42 /* Divide a by b, producing quotient q and remainder r. 43 44 sizeof a is m 45 sizeof b is n 46 sizeof q is m - n 47 sizeof r is n 48 49 The quotient must fit in m - n bytes, i.e., the most significant 50 n digits of a must be less than b, and m must be greater than n. */ 51 52 /* The name of this used to be __div_internal, 53 but that is too long for SYSV. */ 54 55 void 56 __bdiv (a, b, q, r, m, n) 57 unsigned short *a, *b, *q, *r; 58 unsigned long m, n; 59 { 60 unsigned long qhat, rhat; 61 unsigned long acc; 62 unsigned short array[12]; 63 unsigned short *u, *v, *u0, *u1, *u2; 64 unsigned short *v0; 65 int d, qn; 66 int i, j; 67 68 if (m > 16 || n > 8) { 69 #ifdef KERNEL 70 panic("bdiv: out of space"); 71 #else 72 #define EMSG "bdiv: out of space." 73 (void)write(STDERR_FILENO, EMSG, sizeof(EMSG) - 1); 74 abort(); 75 #endif 76 } 77 u = &array[0]; 78 v = &array[m / sizeof(unsigned short)]; 79 m /= sizeof *a; 80 n /= sizeof *b; 81 qn = m - n; 82 83 /* Remove leading zero digits from divisor, and the same number of 84 digits (which must be zero) from dividend. */ 85 86 while (b[big_end (n)] == 0) 87 { 88 r[big_end (n)] = 0; 89 90 a += little_end (2); 91 b += little_end (2); 92 r += little_end (2); 93 m--; 94 n--; 95 96 /* Check for zero divisor. */ 97 if (n == 0) 98 { 99 *r /= n; /* force a divide-by-zero trap */ 100 return; 101 } 102 } 103 104 /* If divisor is a single digit, do short division. */ 105 106 if (n == 1) 107 { 108 acc = a[big_end (m)]; 109 a += little_end (2); 110 for (j = big_end (qn); is_not_lsd (j, qn); j = next_lsd (j)) 111 { 112 acc = (acc << 16) | a[j]; 113 q[j] = acc / *b; 114 acc = acc % *b; 115 } 116 *r = acc; 117 return; 118 } 119 120 /* No such luck, must do long division. Shift divisor and dividend 121 left until the high bit of the divisor is 1. */ 122 123 for (d = 0; d < 16; d++) 124 if (b[big_end (n)] & (1 << (16 - 1 - d))) 125 break; 126 127 bshift (a, d, u, 0, m); 128 bshift (b, d, v, 0, n); 129 130 /* Get pointers to the high dividend and divisor digits. */ 131 132 u0 = u + big_end (m) - big_end (qn); 133 u1 = next_lsd (u0); 134 u2 = next_lsd (u1); 135 u += little_end (2); 136 137 v0 = v + big_end (n); 138 139 /* Main loop: find a quotient digit, multiply it by the divisor, 140 and subtract that from the dividend, shifted over the right amount. */ 141 142 for (j = big_end (qn); is_not_lsd (j, qn); j = next_lsd (j)) 143 { 144 /* Quotient digit initial guess: high 2 dividend digits over high 145 divisor digit. */ 146 147 if (u0[j] == *v0) 148 { 149 qhat = B - 1; 150 rhat = (unsigned long) *v0 + u1[j]; 151 } 152 else 153 { 154 unsigned long numerator = ((unsigned long) u0[j] << 16) | u1[j]; 155 qhat = numerator / *v0; 156 rhat = numerator % *v0; 157 } 158 159 /* Now get the quotient right for high 3 dividend digits over 160 high 2 divisor digits. */ 161 162 while (rhat < B && qhat * *next_lsd (v0) > ((rhat << 16) | u2[j])) 163 { 164 qhat -= 1; 165 rhat += *v0; 166 } 167 168 /* Multiply quotient by divisor, subtract from dividend. */ 169 170 acc = 0; 171 for (i = little_end (n); is_not_msd (i, n); i = next_msd (i)) 172 { 173 acc += (unsigned long) (u + j)[i] - v[i] * qhat; 174 (u + j)[i] = acc & low16; 175 if (acc < B) 176 acc = 0; 177 else 178 acc = (acc >> 16) | -B; 179 } 180 181 q[j] = qhat; 182 183 /* Quotient may have been too high by 1. If dividend went negative, 184 decrement the quotient by 1 and add the divisor back. */ 185 186 if ((signed long) (acc + u0[j]) < 0) 187 { 188 q[j] -= 1; 189 acc = 0; 190 for (i = little_end (n); is_not_msd (i, n); i = next_msd (i)) 191 { 192 acc += (unsigned long) (u + j)[i] + v[i]; 193 (u + j)[i] = acc & low16; 194 acc = acc >> 16; 195 } 196 } 197 } 198 199 /* Now the remainder is what's left of the dividend, shifted right 200 by the amount of the normalizing left shift at the top. */ 201 202 r[big_end (n)] = bshift (u + 1 + little_end (j - 1), 203 16 - d, 204 r + little_end (2), 205 u[little_end (m - 1)] >> d, 206 n - 1); 207 } 208 209 /* Left shift U by K giving W; fill the introduced low-order bits with 210 CARRY_IN. Length of U and W is N. Return carry out. K must be 211 in 0 .. 16. */ 212 213 static int 214 bshift (u, k, w, carry_in, n) 215 unsigned short *u, *w, carry_in; 216 int k, n; 217 { 218 unsigned long acc; 219 int i; 220 221 if (k == 0) 222 { 223 while (n-- > 0) 224 *w++ = *u++; 225 return 0; 226 } 227 228 acc = carry_in; 229 for (i = little_end (n); is_not_msd (i, n); i = next_msd (i)) 230 { 231 acc |= (unsigned long) u[i] << k; 232 w[i] = acc & low16; 233 acc = acc >> 16; 234 } 235 return acc; 236 } 237