-/* original source code from Hackers-Delight
- https://github.com/hcs0/Hackers-Delight
-*/
-/* This divides an n-word dividend by an m-word divisor, giving an
-n-m+1-word quotient and m-word remainder. The bignums are in arrays of
-words. Here a "word" is 32 bits. This routine is designed for a 64-bit
-machine which has a 64/64 division instruction. */
-
-#include <stdio.h>
-#include <stdlib.h> //To define "exit", req'd by XLC.
-#include <stdbool.h>
-#include <stdint.h>
-
-#define max(x, y) ((x) > (y) ? (x) : (y))
-
-int nlz(unsigned x) {
- int n;
-
- if (x == 0) return(32);
- n = 0;
- if (x <= 0x0000FFFF) {n = n +16; x = x <<16;}
- if (x <= 0x00FFFFFF) {n = n + 8; x = x << 8;}
- if (x <= 0x0FFFFFFF) {n = n + 4; x = x << 4;}
- if (x <= 0x3FFFFFFF) {n = n + 2; x = x << 2;}
- if (x <= 0x7FFFFFFF) {n = n + 1;}
- return n;
-}
-
-void dumpit(char *msg, int n, unsigned v[]) {
- int i;
- printf("%s", msg);
- for (i = n-1; i >= 0; i--) printf(" %08x", v[i]);
- printf("\n");
-}
-
-/* q[0], r[0], u[0], and v[0] contain the LEAST significant words.
-(The sequence is in little-endian order).
-
-This is a fairly precise implementation of Knuth's Algorithm D, for a
-binary computer with base b = 2**32. The caller supplies:
- 1. Space q for the quotient, m - n + 1 words (at least one).
- 2. Space r for the remainder (optional), n words.
- 3. The dividend u, m words, m >= 1.
- 4. The divisor v, n words, n >= 2.
-The most significant digit of the divisor, v[n-1], must be nonzero. The
-dividend u may have leading zeros; this just makes the algorithm take
-longer and makes the quotient contain more leading zeros. A value of
-NULL may be given for the address of the remainder to signify that the
-caller does not want the remainder.
- The program does not alter the input parameters u and v.
- The quotient and remainder returned may have leading zeros. The
-function itself returns a value of 0 for success and 1 for invalid
-parameters (e.g., division by 0).
- For now, we must have m >= n. Knuth's Algorithm D also requires
-that the dividend be at least as long as the divisor. (In his terms,
-m >= 0 (unstated). Therefore m+n >= n.) */
-
-int divmnu(unsigned q[], unsigned r[],
- const unsigned u[], const unsigned v[],
- int m, int n) {
-
- const unsigned long long b = 4294967296LL; // Number base (2**32).
- unsigned *un, *vn; // Normalized form of u, v.
- unsigned long long qhat; // Estimated quotient digit.
- unsigned long long rhat; // A remainder.
- unsigned long long p; // Product of two digits.
- long long t, k;
- int s, i, j;
-
- if (m < n || n <= 0 || v[n-1] == 0)
- return 1; // Return if invalid param.
-
- if (n == 1) { // Take care of
- k = 0; // the case of a
- for (j = m - 1; j >= 0; j--) { // single-digit
- q[j] = (k*b + u[j])/v[0]; // divisor here.
- k = (k*b + u[j]) - q[j]*v[0];
- }
- if (r != NULL) r[0] = k;
- return 0;
- }
-
- /* Normalize by shifting v left just enough so that its high-order
- bit is on, and shift u left the same amount. We may have to append a
- high-order digit on the dividend; we do that unconditionally. */
-
- s = nlz(v[n-1]); // 0 <= s <= 31.
- vn = (unsigned *)alloca(4*n);
- for (i = n - 1; i > 0; i--)
- vn[i] = (v[i] << s) | ((unsigned long long)v[i-1] >> (32-s));
- vn[0] = v[0] << s;
-
- un = (unsigned *)alloca(4*(m + 1));
- un[m] = (unsigned long long)u[m-1] >> (32-s);
- for (i = m - 1; i > 0; i--)
- un[i] = (u[i] << s) | ((unsigned long long)u[i-1] >> (32-s));
- un[0] = u[0] << s;
-
- for (j = m - n; j >= 0; j--) { // Main loop.
- // Compute estimate qhat of q[j].
- qhat = (un[j+n]*b + un[j+n-1])/vn[n-1];
- rhat = (un[j+n]*b + un[j+n-1]) - qhat*vn[n-1];
-again:
- if (qhat >= b || qhat*vn[n-2] > b*rhat + un[j+n-2])
- { qhat = qhat - 1;
- rhat = rhat + vn[n-1];
- if (rhat < b) goto again;
- }
-
-#define SUB_MUL_BORROW
-#ifdef ORIGINAL
- // Multiply and subtract.
- k = 0;
- for (i = 0; i < n; i++) {
- p = qhat*vn[i];
- t = un[i+j] - k - (p & 0xFFFFFFFFLL);
- un[i+j] = t;
- k = (p >> 32) - (t >> 32);
- }
- t = un[j+n] - k;
- un[j+n] = t;
- bool need_fixup = t < 0;
-#elif defined(SUB_MUL_BORROW)
- (void)p; // shut up unused variable warning
-
- // Multiply and subtract.
- uint32_t borrow = 0;
- uint32_t phi[2000]; // plenty space
- uint32_t plo[2000]; // plenty space
- // first, perform mul-and-sub and store in split hi-lo
- // this shows the vectorised sv.msubx which stores 128-bit in
- // two 64-bit registers
- for(int i = 0; i <= n; i++) {
- uint32_t vn_i = i < n ? vn[i] : 0;
- uint64_t value = un[i + j] - (uint64_t)qhat * vn_i;
- plo[i] = value & 0xffffffffLL;
- phi[i] = value >> 32;
- }
- // second, reconstruct the 64-bit result, subtract borrow,
- // store top-half (-ve) in new borrow and store low-half as answer
- // this is the new (odd) instruction
- for(int i = 0; i <= n; i++) {
- uint64_t value = (((uint64_t)phi[i]<<32) | plo[i]) - borrow;
- borrow = -(uint32_t)(value >> 32);
- un[i + j] = (uint32_t)value;
- }
- bool need_fixup = borrow != 0;
-#elif defined(MUL_RSUB_CARRY)
- (void)p; // shut up unused variable warning
-
- // Multiply and subtract.
- uint32_t carry = 1;
- for(int i = 0; i <= n; i++) {
- uint32_t vn_i = i < n ? vn[i] : 0;
- uint64_t result = un[i + j] + ~((uint64_t)qhat * vn_i) + carry;
- uint32_t result_high = result >> 32;
- if(carry <= 1)
- result_high++;
- carry = result_high;
- un[i + j] = (uint32_t)result;
- }
- bool need_fixup = carry != 1;
-#else
-#error need to define one of ORIGINAL, SUB_MUL_BORROW, or MUL_RSUB_CARRY
-#endif
-
- q[j] = qhat; // Store quotient digit.
- if (need_fixup) { // If we subtracted too
- q[j] = q[j] - 1; // much, add back.
- k = 0;
- for (i = 0; i < n; i++) {
- t = (unsigned long long)un[i+j] + vn[i] + k;
- un[i+j] = t;
- k = t >> 32;
- }
- un[j+n] = un[j+n] + k;
- }
- } // End j.
- // If the caller wants the remainder, unnormalize
- // it and pass it back.
- if (r != NULL) {
- for (i = 0; i < n-1; i++)
- r[i] = (un[i] >> s) | ((unsigned long long)un[i+1] << (32-s));
- r[n-1] = un[n-1] >> s;
- }
- return 0;
-}
-
-int errors;
-
-void check(unsigned q[], unsigned r[],
- unsigned u[], unsigned v[],
- int m, int n,
- unsigned cq[], unsigned cr[]) {
- int i, szq;
-
- szq = max(m - n + 1, 1);
- for (i = 0; i < szq; i++) {
- if (q[i] != cq[i]) {
- errors = errors + 1;
- dumpit("Error, dividend u =", m, u);
- dumpit(" divisor v =", n, v);
- dumpit("For quotient, got:", m-n+1, q);
- dumpit(" Should get:", m-n+1, cq);
- return;
- }
- }
- for (i = 0; i < n; i++) {
- if (r[i] != cr[i]) {
- errors = errors + 1;
- dumpit("Error, dividend u =", m, u);
- dumpit(" divisor v =", n, v);
- dumpit("For remainder, got:", n, r);
- dumpit(" Should get:", n, cr);
- return;
- }
- }
- return;
-}
-
-int main() {
- static unsigned test[] = {
- // m, n, u..., v..., cq..., cr....
- 1, 1, 3, 0, 1, 1, // Error, divide by 0.
- 1, 2, 7, 1,3, 0, 7,0, // Error, n > m.
- 2, 2, 0,0, 1,0, 0, 0,0, // Error, incorrect remainder cr.
- 1, 1, 3, 2, 1, 1,
- 1, 1, 3, 3, 1, 0,
- 1, 1, 3, 4, 0, 3,
- 1, 1, 0, 0xffffffff, 0, 0,
- 1, 1, 0xffffffff, 1, 0xffffffff, 0,
- 1, 1, 0xffffffff, 0xffffffff, 1, 0,
- 1, 1, 0xffffffff, 3, 0x55555555, 0,
- 2, 1, 0xffffffff,0xffffffff, 1, 0xffffffff,0xffffffff, 0,
- 2, 1, 0xffffffff,0xffffffff, 0xffffffff, 1,1, 0,
- 2, 1, 0xffffffff,0xfffffffe, 0xffffffff, 0xffffffff,0, 0xfffffffe,
- 2, 1, 0x00005678,0x00001234, 0x00009abc, 0x1e1dba76,0, 0x6bd0,
- 2, 2, 0,0, 0,1, 0, 0,0,
- 2, 2, 0,7, 0,3, 2, 0,1,
- 2, 2, 5,7, 0,3, 2, 5,1,
- 2, 2, 0,6, 0,2, 3, 0,0,
- 1, 1, 0x80000000, 0x40000001, 0x00000001, 0x3fffffff,
- 2, 1, 0x00000000,0x80000000, 0x40000001, 0xfffffff8,0x00000001, 0x00000008,
- 2, 2, 0x00000000,0x80000000, 0x00000001,0x40000000, 0x00000001, 0xffffffff,0x3fffffff,
- 2, 2, 0x0000789a,0x0000bcde, 0x0000789a,0x0000bcde, 1, 0,0,
- 2, 2, 0x0000789b,0x0000bcde, 0x0000789a,0x0000bcde, 1, 1,0,
- 2, 2, 0x00007899,0x0000bcde, 0x0000789a,0x0000bcde, 0, 0x00007899,0x0000bcde,
- 2, 2, 0x0000ffff,0x0000ffff, 0x0000ffff,0x0000ffff, 1, 0,0,
- 2, 2, 0x0000ffff,0x0000ffff, 0x00000000,0x00000001, 0x0000ffff, 0x0000ffff,0,
- 3, 2, 0x000089ab,0x00004567,0x00000123, 0x00000000,0x00000001, 0x00004567,0x00000123, 0x000089ab,0,
- 3, 2, 0x00000000,0x0000fffe,0x00008000, 0x0000ffff,0x00008000, 0xffffffff,0x00000000, 0x0000ffff,0x00007fff, // Shows that first qhat can = b + 1.
- 3, 3, 0x00000003,0x00000000,0x80000000, 0x00000001,0x00000000,0x20000000, 0x00000003, 0,0,0x20000000, // Adding back step req'd.
- 3, 3, 0x00000003,0x00000000,0x00008000, 0x00000001,0x00000000,0x00002000, 0x00000003, 0,0,0x00002000, // Adding back step req'd.
- 4, 3, 0,0,0x00008000,0x00007fff, 1,0,0x00008000, 0xfffe0000,0, 0x00020000,0xffffffff,0x00007fff, // Add back req'd.
- 4, 3, 0,0x0000fffe,0,0x00008000, 0x0000ffff,0,0x00008000, 0xffffffff,0, 0x0000ffff,0xffffffff,0x00007fff, // Shows that mult-sub quantity cannot be treated as signed.
- 4, 3, 0,0xfffffffe,0,0x80000000, 0x0000ffff,0,0x80000000, 0x00000000,1, 0x00000000,0xfffeffff,0x00000000, // Shows that mult-sub quantity cannot be treated as signed.
- 4, 3, 0,0xfffffffe,0,0x80000000, 0xffffffff,0,0x80000000, 0xffffffff,0, 0xffffffff,0xffffffff,0x7fffffff, // Shows that mult-sub quantity cannot be treated as signed.
- };
- int i, n, m, ncases, f;
- unsigned q[10], r[10];
- unsigned *u, *v, *cq, *cr;
-
- printf("divmnu:\n");
- i = 0;
- ncases = 0;
- while (i < sizeof(test)/4) {
- m = test[i];
- n = test[i+1];
- u = &test[i+2];
- v = &test[i+2+m];
- cq = &test[i+2+m+n];
- cr = &test[i+2+m+n+max(m-n+1, 1)];
-
- f = divmnu(q, r, u, v, m, n);
- if (f) {
- dumpit("Error return code for dividend u =", m, u);
- dumpit(" divisor v =", n, v);
- errors = errors + 1;
- }
- else
- check(q, r, u, v, m, n, cq, cr);
- i = i + 2 + m + n + max(m-n+1, 1) + n;
- ncases = ncases + 1;
- }
-
- printf("%d errors out of %d cases; there should be 3.\n", errors, ncases);
- return 0;
-}
+// moved to ../biginteger/divmnu.c
projects / libreriscv.git / blobdiff