Commit | Line | Data |
---|---|---|
90b1a5c8 KB |
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) %G%"; | |
10 | #endif /* LIBC_SCCS and not lint */ | |
11 | ||
ce3a66ab KM |
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) | |
c974bcd2 KM |
62 | { |
63 | *r /= n; /* force a divide-by-zero trap */ | |
64 | return; | |
65 | } | |
ce3a66ab KM |
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 | { | |
c974bcd2 KM |
187 | while (n-- > 0) |
188 | *w++ = *u++; | |
ce3a66ab KM |
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 | } |