1 /* original source code from Hackers-Delight
2 https://github.com/hcs0/Hackers-Delight
4 /* This divides an n-word dividend by an m-word divisor, giving an
5 n-m+1-word quotient and m-word remainder. The bignums are in arrays of
6 words. Here a "word" is 32 bits. This routine is designed for a 64-bit
7 machine which has a 64/64 division instruction. */
10 #include <stdlib.h> //To define "exit", req'd by XLC.
14 #define max(x, y) ((x) > (y) ? (x) : (y))
19 if (x
== 0) return(32);
21 if (x
<= 0x0000FFFF) {n
= n
+16; x
= x
<<16;}
22 if (x
<= 0x00FFFFFF) {n
= n
+ 8; x
= x
<< 8;}
23 if (x
<= 0x0FFFFFFF) {n
= n
+ 4; x
= x
<< 4;}
24 if (x
<= 0x3FFFFFFF) {n
= n
+ 2; x
= x
<< 2;}
25 if (x
<= 0x7FFFFFFF) {n
= n
+ 1;}
29 void dumpit(char *msg
, int n
, unsigned v
[]) {
32 for (i
= n
-1; i
>= 0; i
--) printf(" %08x", v
[i
]);
36 /* q[0], r[0], u[0], and v[0] contain the LEAST significant words.
37 (The sequence is in little-endian order).
39 This is a fairly precise implementation of Knuth's Algorithm D, for a
40 binary computer with base b = 2**32. The caller supplies:
41 1. Space q for the quotient, m - n + 1 words (at least one).
42 2. Space r for the remainder (optional), n words.
43 3. The dividend u, m words, m >= 1.
44 4. The divisor v, n words, n >= 2.
45 The most significant digit of the divisor, v[n-1], must be nonzero. The
46 dividend u may have leading zeros; this just makes the algorithm take
47 longer and makes the quotient contain more leading zeros. A value of
48 NULL may be given for the address of the remainder to signify that the
49 caller does not want the remainder.
50 The program does not alter the input parameters u and v.
51 The quotient and remainder returned may have leading zeros. The
52 function itself returns a value of 0 for success and 1 for invalid
53 parameters (e.g., division by 0).
54 For now, we must have m >= n. Knuth's Algorithm D also requires
55 that the dividend be at least as long as the divisor. (In his terms,
56 m >= 0 (unstated). Therefore m+n >= n.) */
58 int divmnu(unsigned q
[], unsigned r
[],
59 const unsigned u
[], const unsigned v
[],
62 const unsigned long long b
= 4294967296LL; // Number base (2**32).
63 unsigned *un
, *vn
; // Normalized form of u, v.
64 unsigned long long qhat
; // Estimated quotient digit.
65 unsigned long long rhat
; // A remainder.
66 unsigned long long p
; // Product of two digits.
70 if (m
< n
|| n
<= 0 || v
[n
-1] == 0)
71 return 1; // Return if invalid param.
73 if (n
== 1) { // Take care of
74 k
= 0; // the case of a
75 for (j
= m
- 1; j
>= 0; j
--) { // single-digit
76 q
[j
] = (k
*b
+ u
[j
])/v
[0]; // divisor here.
77 k
= (k
*b
+ u
[j
]) - q
[j
]*v
[0];
79 if (r
!= NULL
) r
[0] = k
;
83 /* Normalize by shifting v left just enough so that its high-order
84 bit is on, and shift u left the same amount. We may have to append a
85 high-order digit on the dividend; we do that unconditionally. */
87 s
= nlz(v
[n
-1]); // 0 <= s <= 31.
88 vn
= (unsigned *)alloca(4*n
);
89 for (i
= n
- 1; i
> 0; i
--)
90 vn
[i
] = (v
[i
] << s
) | ((unsigned long long)v
[i
-1] >> (32-s
));
93 un
= (unsigned *)alloca(4*(m
+ 1));
94 un
[m
] = (unsigned long long)u
[m
-1] >> (32-s
);
95 for (i
= m
- 1; i
> 0; i
--)
96 un
[i
] = (u
[i
] << s
) | ((unsigned long long)u
[i
-1] >> (32-s
));
99 for (j
= m
- n
; j
>= 0; j
--) { // Main loop.
100 // Compute estimate qhat of q[j].
101 qhat
= (un
[j
+n
]*b
+ un
[j
+n
-1])/vn
[n
-1];
102 rhat
= (un
[j
+n
]*b
+ un
[j
+n
-1]) - qhat
*vn
[n
-1];
104 if (qhat
>= b
|| qhat
*vn
[n
-2] > b
*rhat
+ un
[j
+n
-2])
106 rhat
= rhat
+ vn
[n
-1];
107 if (rhat
< b
) goto again
;
111 // Multiply and subtract.
113 for (i
= 0; i
< n
; i
++) {
115 t
= un
[i
+j
] - k
- (p
& 0xFFFFFFFFLL
);
117 k
= (p
>> 32) - (t
>> 32);
121 bool need_fixup
= t
< 0;
122 #elif defined(SUB_MUL_BORROW)
123 (void)p
; // shut up unused variable warning
125 // Multiply and subtract.
127 for(int i
= 0; i
<= n
; i
++) {
128 uint32_t vn_i
= i
< n
? vn
[i
] : 0;
129 uint64_t value
= un
[i
+ j
] - (uint64_t)qhat
* vn_i
- borrow
;
130 borrow
= -(uint32_t)(value
>> 32);
131 un
[i
+ j
] = (uint32_t)value
;
133 bool need_fixup
= borrow
!= 0;
134 #elif defined(MUL_RSUB_CARRY)
135 (void)p
; // shut up unused variable warning
137 // Multiply and subtract.
139 for(int i
= 0; i
<= n
; i
++) {
140 uint32_t vn_i
= i
< n
? vn
[i
] : 0;
141 uint64_t result
= un
[i
+ j
] + ~((uint64_t)qhat
* vn_i
) + carry
;
142 uint32_t result_high
= result
>> 32;
146 un
[i
+ j
] = (uint32_t)result
;
148 bool need_fixup
= carry
!= 1;
149 #elif defined(SUB_MUL_BORROW_2_STAGE)
150 (void)p
; // shut up unused variable warning
152 // Multiply and subtract.
154 uint32_t phi
[2000]; // plenty space
155 uint32_t plo
[2000]; // plenty space
156 // first, perform mul-and-sub and store in split hi-lo
157 // this shows the vectorised sv.msubx which stores 128-bit in
158 // two 64-bit registers
159 for(int i
= 0; i
<= n
; i
++) {
160 uint32_t vn_i
= i
< n
? vn
[i
] : 0;
161 uint64_t value
= un
[i
+ j
] - (uint64_t)qhat
* vn_i
;
162 plo
[i
] = value
& 0xffffffffLL
;
163 phi
[i
] = value
>> 32;
165 // second, reconstruct the 64-bit result, subtract borrow,
166 // store top-half (-ve) in new borrow and store low-half as answer
167 // this is the new (odd) instruction
168 for(int i
= 0; i
<= n
; i
++) {
169 uint64_t value
= (((uint64_t)phi
[i
]<<32) | plo
[i
]) - borrow
;
170 borrow
= ~(value
>> 32)+1; // -(uint32_t)(value >> 32);
171 un
[i
+ j
] = (uint32_t)value
;
173 bool need_fixup
= borrow
!= 0;
174 #elif defined(MUL_RSUB_CARRY_2_STAGE)
175 (void)p
; // shut up unused variable warning
177 // Multiply and subtract.
179 uint32_t phi
[2000]; // plenty space
180 uint32_t plo
[2000]; // plenty space
181 for(int i
= 0; i
<= n
; i
++) {
182 uint32_t vn_i
= i
< n
? vn
[i
] : 0;
183 uint64_t value
= un
[i
+ j
] + ~((uint64_t)qhat
* vn_i
);
184 plo
[i
] = value
& 0xffffffffLL
;
185 phi
[i
] = value
>> 32;
187 for(int i
= 0; i
<= n
; i
++) {
188 uint64_t result
= (((uint64_t)phi
[i
]<<32) | plo
[i
]) + carry
;
189 uint32_t result_high
= result
>> 32;
193 un
[i
+ j
] = (uint32_t)result
;
195 bool need_fixup
= carry
!= 1;
196 #elif defined(MUL_RSUB_CARRY_2_STAGE1)
197 (void)p
; // shut up unused variable warning
199 // Multiply and subtract.
201 uint32_t phi
[2000]; // plenty space
202 uint32_t plo
[2000]; // plenty space
203 // same mul-and-sub as SUB_MUL_BORROW but not the same
204 // mul-and-sub-minus-one as MUL_RSUB_CARRY
205 for(int i
= 0; i
<= n
; i
++) {
206 uint32_t vn_i
= i
< n
? vn
[i
] : 0;
207 uint64_t value
= un
[i
+ j
] - ((uint64_t)qhat
* vn_i
);
208 plo
[i
] = value
& 0xffffffffLL
;
209 phi
[i
] = value
>> 32;
211 // compensate for the +1 that was added by mul-and-sub by subtracting
213 for(int i
= 0; i
<= n
; i
++) {
214 uint64_t result
= (((uint64_t)phi
[i
]<<32) | plo
[i
]) + carry
+
215 ~(0); // a way to express "-1"
216 uint32_t result_high
= result
>> 32;
220 un
[i
+ j
] = (uint32_t)result
;
222 bool need_fixup
= carry
!= 1;
224 #error need to choose one of the algorithm options; e.g. -DORIGINAL
227 q
[j
] = qhat
; // Store quotient digit.
228 if (need_fixup
) { // If we subtracted too
229 q
[j
] = q
[j
] - 1; // much, add back.
231 for (i
= 0; i
< n
; i
++) {
232 t
= (unsigned long long)un
[i
+j
] + vn
[i
] + k
;
236 un
[j
+n
] = un
[j
+n
] + k
;
239 // If the caller wants the remainder, unnormalize
240 // it and pass it back.
242 for (i
= 0; i
< n
-1; i
++)
243 r
[i
] = (un
[i
] >> s
) | ((unsigned long long)un
[i
+1] << (32-s
));
244 r
[n
-1] = un
[n
-1] >> s
;
251 void check(unsigned q
[], unsigned r
[],
252 unsigned u
[], unsigned v
[],
254 unsigned cq
[], unsigned cr
[]) {
257 szq
= max(m
- n
+ 1, 1);
258 for (i
= 0; i
< szq
; i
++) {
261 dumpit("Error, dividend u =", m
, u
);
262 dumpit(" divisor v =", n
, v
);
263 dumpit("For quotient, got:", m
-n
+1, q
);
264 dumpit(" Should get:", m
-n
+1, cq
);
268 for (i
= 0; i
< n
; i
++) {
271 dumpit("Error, dividend u =", m
, u
);
272 dumpit(" divisor v =", n
, v
);
273 dumpit("For remainder, got:", n
, r
);
274 dumpit(" Should get:", n
, cr
);
282 static unsigned test
[] = {
283 // m, n, u..., v..., cq..., cr....
284 1, 1, 3, 0, 1, 1, // Error, divide by 0.
285 1, 2, 7, 1,3, 0, 7,0, // Error, n > m.
286 2, 2, 0,0, 1,0, 0, 0,0, // Error, incorrect remainder cr.
290 1, 1, 0, 0xffffffff, 0, 0,
291 1, 1, 0xffffffff, 1, 0xffffffff, 0,
292 1, 1, 0xffffffff, 0xffffffff, 1, 0,
293 1, 1, 0xffffffff, 3, 0x55555555, 0,
294 2, 1, 0xffffffff,0xffffffff, 1, 0xffffffff,0xffffffff, 0,
295 2, 1, 0xffffffff,0xffffffff, 0xffffffff, 1,1, 0,
296 2, 1, 0xffffffff,0xfffffffe, 0xffffffff, 0xffffffff,0, 0xfffffffe,
297 2, 1, 0x00005678,0x00001234, 0x00009abc, 0x1e1dba76,0, 0x6bd0,
298 2, 2, 0,0, 0,1, 0, 0,0,
299 2, 2, 0,7, 0,3, 2, 0,1,
300 2, 2, 5,7, 0,3, 2, 5,1,
301 2, 2, 0,6, 0,2, 3, 0,0,
302 1, 1, 0x80000000, 0x40000001, 0x00000001, 0x3fffffff,
303 2, 1, 0x00000000,0x80000000, 0x40000001, 0xfffffff8,0x00000001, 0x00000008,
304 2, 2, 0x00000000,0x80000000, 0x00000001,0x40000000, 0x00000001, 0xffffffff,0x3fffffff,
305 2, 2, 0x0000789a,0x0000bcde, 0x0000789a,0x0000bcde, 1, 0,0,
306 2, 2, 0x0000789b,0x0000bcde, 0x0000789a,0x0000bcde, 1, 1,0,
307 2, 2, 0x00007899,0x0000bcde, 0x0000789a,0x0000bcde, 0, 0x00007899,0x0000bcde,
308 2, 2, 0x0000ffff,0x0000ffff, 0x0000ffff,0x0000ffff, 1, 0,0,
309 2, 2, 0x0000ffff,0x0000ffff, 0x00000000,0x00000001, 0x0000ffff, 0x0000ffff,0,
310 3, 2, 0x000089ab,0x00004567,0x00000123, 0x00000000,0x00000001, 0x00004567,0x00000123, 0x000089ab,0,
311 3, 2, 0x00000000,0x0000fffe,0x00008000, 0x0000ffff,0x00008000, 0xffffffff,0x00000000, 0x0000ffff,0x00007fff, // Shows that first qhat can = b + 1.
312 3, 3, 0x00000003,0x00000000,0x80000000, 0x00000001,0x00000000,0x20000000, 0x00000003, 0,0,0x20000000, // Adding back step req'd.
313 3, 3, 0x00000003,0x00000000,0x00008000, 0x00000001,0x00000000,0x00002000, 0x00000003, 0,0,0x00002000, // Adding back step req'd.
314 4, 3, 0,0,0x00008000,0x00007fff, 1,0,0x00008000, 0xfffe0000,0, 0x00020000,0xffffffff,0x00007fff, // Add back req'd.
315 4, 3, 0,0x0000fffe,0,0x00008000, 0x0000ffff,0,0x00008000, 0xffffffff,0, 0x0000ffff,0xffffffff,0x00007fff, // Shows that mult-sub quantity cannot be treated as signed.
316 4, 3, 0,0xfffffffe,0,0x80000000, 0x0000ffff,0,0x80000000, 0x00000000,1, 0x00000000,0xfffeffff,0x00000000, // Shows that mult-sub quantity cannot be treated as signed.
317 4, 3, 0,0xfffffffe,0,0x80000000, 0xffffffff,0,0x80000000, 0xffffffff,0, 0xffffffff,0xffffffff,0x7fffffff, // Shows that mult-sub quantity cannot be treated as signed.
319 int i
, n
, m
, ncases
, f
;
320 unsigned q
[10], r
[10];
321 unsigned *u
, *v
, *cq
, *cr
;
326 while (i
< sizeof(test
)/4) {
332 cr
= &test
[i
+2+m
+n
+max(m
-n
+1, 1)];
334 f
= divmnu(q
, r
, u
, v
, m
, n
);
336 dumpit("Error return code for dividend u =", m
, u
);
337 dumpit(" divisor v =", n
, v
);
341 check(q
, r
, u
, v
, m
, n
, cq
, cr
);
342 i
= i
+ 2 + m
+ n
+ max(m
-n
+1, 1) + n
;
346 printf("%d errors out of %d cases; there should be 3.\n", errors
, ncases
);