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
;
110 #define SUB_MUL_BORROW
112 // Multiply and subtract.
114 for (i
= 0; i
< n
; i
++) {
116 t
= un
[i
+j
] - k
- (p
& 0xFFFFFFFFLL
);
118 k
= (p
>> 32) - (t
>> 32);
122 bool need_fixup
= t
< 0;
123 #elif defined(SUB_MUL_BORROW)
124 (void)p
; // shut up unused variable warning
126 // Multiply and subtract.
128 uint32_t phi
[2000]; // plenty space
129 uint32_t plo
[2000]; // plenty space
130 // first, perform mul-and-sub and store in split hi-lo
131 // this shows the vectorised sv.msubx which stores 128-bit in
132 // two 64-bit registers
133 for(int i
= 0; i
<= n
; i
++) {
134 uint32_t vn_i
= i
< n
? vn
[i
] : 0;
135 uint64_t value
= un
[i
+ j
] - (uint64_t)qhat
* vn_i
;
136 plo
[i
] = value
& 0xffffffffLL
;
137 phi
[i
] = value
>> 32;
139 // second, reconstruct the 64-bit result, subtract borrow,
140 // store top-half (-ve) in new borrow and store low-half as answer
141 // this is the new (odd) instruction
142 for(int i
= 0; i
<= n
; i
++) {
143 uint64_t value
= (((uint64_t)phi
[i
]<<32) | plo
[i
]) - borrow
;
144 borrow
= -(uint32_t)(value
>> 32);
145 un
[i
+ j
] = (uint32_t)value
;
147 bool need_fixup
= borrow
!= 0;
148 #elif defined(MUL_RSUB_CARRY)
149 (void)p
; // shut up unused variable warning
151 // Multiply and subtract.
153 for(int i
= 0; i
<= n
; i
++) {
154 uint32_t vn_i
= i
< n
? vn
[i
] : 0;
155 uint64_t result
= un
[i
+ j
] + ~((uint64_t)qhat
* vn_i
) + carry
;
156 uint32_t result_high
= result
>> 32;
160 un
[i
+ j
] = (uint32_t)result
;
162 bool need_fixup
= carry
!= 1;
164 #error need to define one of ORIGINAL, SUB_MUL_BORROW, or MUL_RSUB_CARRY
167 q
[j
] = qhat
; // Store quotient digit.
168 if (need_fixup
) { // If we subtracted too
169 q
[j
] = q
[j
] - 1; // much, add back.
171 for (i
= 0; i
< n
; i
++) {
172 t
= (unsigned long long)un
[i
+j
] + vn
[i
] + k
;
176 un
[j
+n
] = un
[j
+n
] + k
;
179 // If the caller wants the remainder, unnormalize
180 // it and pass it back.
182 for (i
= 0; i
< n
-1; i
++)
183 r
[i
] = (un
[i
] >> s
) | ((unsigned long long)un
[i
+1] << (32-s
));
184 r
[n
-1] = un
[n
-1] >> s
;
191 void check(unsigned q
[], unsigned r
[],
192 unsigned u
[], unsigned v
[],
194 unsigned cq
[], unsigned cr
[]) {
197 szq
= max(m
- n
+ 1, 1);
198 for (i
= 0; i
< szq
; i
++) {
201 dumpit("Error, dividend u =", m
, u
);
202 dumpit(" divisor v =", n
, v
);
203 dumpit("For quotient, got:", m
-n
+1, q
);
204 dumpit(" Should get:", m
-n
+1, cq
);
208 for (i
= 0; i
< n
; i
++) {
211 dumpit("Error, dividend u =", m
, u
);
212 dumpit(" divisor v =", n
, v
);
213 dumpit("For remainder, got:", n
, r
);
214 dumpit(" Should get:", n
, cr
);
222 static unsigned test
[] = {
223 // m, n, u..., v..., cq..., cr....
224 1, 1, 3, 0, 1, 1, // Error, divide by 0.
225 1, 2, 7, 1,3, 0, 7,0, // Error, n > m.
226 2, 2, 0,0, 1,0, 0, 0,0, // Error, incorrect remainder cr.
230 1, 1, 0, 0xffffffff, 0, 0,
231 1, 1, 0xffffffff, 1, 0xffffffff, 0,
232 1, 1, 0xffffffff, 0xffffffff, 1, 0,
233 1, 1, 0xffffffff, 3, 0x55555555, 0,
234 2, 1, 0xffffffff,0xffffffff, 1, 0xffffffff,0xffffffff, 0,
235 2, 1, 0xffffffff,0xffffffff, 0xffffffff, 1,1, 0,
236 2, 1, 0xffffffff,0xfffffffe, 0xffffffff, 0xffffffff,0, 0xfffffffe,
237 2, 1, 0x00005678,0x00001234, 0x00009abc, 0x1e1dba76,0, 0x6bd0,
238 2, 2, 0,0, 0,1, 0, 0,0,
239 2, 2, 0,7, 0,3, 2, 0,1,
240 2, 2, 5,7, 0,3, 2, 5,1,
241 2, 2, 0,6, 0,2, 3, 0,0,
242 1, 1, 0x80000000, 0x40000001, 0x00000001, 0x3fffffff,
243 2, 1, 0x00000000,0x80000000, 0x40000001, 0xfffffff8,0x00000001, 0x00000008,
244 2, 2, 0x00000000,0x80000000, 0x00000001,0x40000000, 0x00000001, 0xffffffff,0x3fffffff,
245 2, 2, 0x0000789a,0x0000bcde, 0x0000789a,0x0000bcde, 1, 0,0,
246 2, 2, 0x0000789b,0x0000bcde, 0x0000789a,0x0000bcde, 1, 1,0,
247 2, 2, 0x00007899,0x0000bcde, 0x0000789a,0x0000bcde, 0, 0x00007899,0x0000bcde,
248 2, 2, 0x0000ffff,0x0000ffff, 0x0000ffff,0x0000ffff, 1, 0,0,
249 2, 2, 0x0000ffff,0x0000ffff, 0x00000000,0x00000001, 0x0000ffff, 0x0000ffff,0,
250 3, 2, 0x000089ab,0x00004567,0x00000123, 0x00000000,0x00000001, 0x00004567,0x00000123, 0x000089ab,0,
251 3, 2, 0x00000000,0x0000fffe,0x00008000, 0x0000ffff,0x00008000, 0xffffffff,0x00000000, 0x0000ffff,0x00007fff, // Shows that first qhat can = b + 1.
252 3, 3, 0x00000003,0x00000000,0x80000000, 0x00000001,0x00000000,0x20000000, 0x00000003, 0,0,0x20000000, // Adding back step req'd.
253 3, 3, 0x00000003,0x00000000,0x00008000, 0x00000001,0x00000000,0x00002000, 0x00000003, 0,0,0x00002000, // Adding back step req'd.
254 4, 3, 0,0,0x00008000,0x00007fff, 1,0,0x00008000, 0xfffe0000,0, 0x00020000,0xffffffff,0x00007fff, // Add back req'd.
255 4, 3, 0,0x0000fffe,0,0x00008000, 0x0000ffff,0,0x00008000, 0xffffffff,0, 0x0000ffff,0xffffffff,0x00007fff, // Shows that mult-sub quantity cannot be treated as signed.
256 4, 3, 0,0xfffffffe,0,0x80000000, 0x0000ffff,0,0x80000000, 0x00000000,1, 0x00000000,0xfffeffff,0x00000000, // Shows that mult-sub quantity cannot be treated as signed.
257 4, 3, 0,0xfffffffe,0,0x80000000, 0xffffffff,0,0x80000000, 0xffffffff,0, 0xffffffff,0xffffffff,0x7fffffff, // Shows that mult-sub quantity cannot be treated as signed.
259 int i
, n
, m
, ncases
, f
;
260 unsigned q
[10], r
[10];
261 unsigned *u
, *v
, *cq
, *cr
;
266 while (i
< sizeof(test
)/4) {
272 cr
= &test
[i
+2+m
+n
+max(m
-n
+1, 1)];
274 f
= divmnu(q
, r
, u
, v
, m
, n
);
276 dumpit("Error return code for dividend u =", m
, u
);
277 dumpit(" divisor v =", n
, v
);
281 check(q
, r
, u
, v
, m
, n
, cq
, cr
);
282 i
= i
+ 2 + m
+ n
+ max(m
-n
+1, 1) + n
;
286 printf("%d errors out of %d cases; there should be 3.\n", errors
, ncases
);