openpower/sv/bitmanip/divmnu64.c

   1 /* original source code from Hackers-Delight
   2    https://github.com/hcs0/Hackers-Delight
   3 */
   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. */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>     //To define "exit", req'd by XLC.
  11
  12 #define max(x, y) ((x) > (y) ? (x) : (y))
  13
  14 int nlz(unsigned x) {
  15    int n;
  16
  17    if (x == 0) return(32);
  18    n = 0;
  19    if (x <= 0x0000FFFF) {n = n +16; x = x <<16;}
  20    if (x <= 0x00FFFFFF) {n = n + 8; x = x << 8;}
  21    if (x <= 0x0FFFFFFF) {n = n + 4; x = x << 4;}
  22    if (x <= 0x3FFFFFFF) {n = n + 2; x = x << 2;}
  23    if (x <= 0x7FFFFFFF) {n = n + 1;}
  24    return n;
  25 }
  26
  27 void dumpit(char *msg, int n, unsigned v[]) {
  28    int i;
  29    printf(msg);
  30    for (i = n-1; i >= 0; i--) printf(" %08x", v[i]);
  31    printf("\n");
  32 }
  33
  34 /* q[0], r[0], u[0], and v[0] contain the LEAST significant words.
  35 (The sequence is in little-endian order).
  36
  37 This is a fairly precise implementation of Knuth's Algorithm D, for a
  38 binary computer with base b = 2**32. The caller supplies:
  39    1. Space q for the quotient, m - n + 1 words (at least one).
  40    2. Space r for the remainder (optional), n words.
  41    3. The dividend u, m words, m >= 1.
  42    4. The divisor v, n words, n >= 2.
  43 The most significant digit of the divisor, v[n-1], must be nonzero.  The
  44 dividend u may have leading zeros; this just makes the algorithm take
  45 longer and makes the quotient contain more leading zeros.  A value of
  46 NULL may be given for the address of the remainder to signify that the
  47 caller does not want the remainder.
  48    The program does not alter the input parameters u and v.
  49    The quotient and remainder returned may have leading zeros.  The
  50 function itself returns a value of 0 for success and 1 for invalid
  51 parameters (e.g., division by 0).
  52    For now, we must have m >= n.  Knuth's Algorithm D also requires
  53 that the dividend be at least as long as the divisor.  (In his terms,
  54 m >= 0 (unstated).  Therefore m+n >= n.) */
  55
  56 int divmnu(unsigned q[], unsigned r[],
  57      const unsigned u[], const unsigned v[],
  58      int m, int n) {
  59
  60    const unsigned long long b = 4294967296LL; // Number base (2**32).
  61    unsigned *un, *vn;                         // Normalized form of u, v.
  62    unsigned long long qhat;                   // Estimated quotient digit.
  63    unsigned long long rhat;                   // A remainder.
  64    unsigned long long p;                      // Product of two digits.
  65    long long t, k;
  66    int s, i, j;
  67
  68    if (m < n || n <= 0 || v[n-1] == 0)
  69       return 1;                         // Return if invalid param.
  70
  71    if (n == 1) {                        // Take care of
  72       k = 0;                            // the case of a
  73       for (j = m - 1; j >= 0; j--) {    // single-digit
  74          q[j] = (k*b + u[j])/v[0];      // divisor here.
  75          k = (k*b + u[j]) - q[j]*v[0];
  76       }
  77       if (r != NULL) r[0] = k;
  78       return 0;
  79    }
  80
  81    /* Normalize by shifting v left just enough so that its high-order
  82    bit is on, and shift u left the same amount. We may have to append a
  83    high-order digit on the dividend; we do that unconditionally. */
  84
  85    s = nlz(v[n-1]);             // 0 <= s <= 31.
  86    vn = (unsigned *)alloca(4*n);
  87    for (i = n - 1; i > 0; i--)
  88       vn[i] = (v[i] << s) | ((unsigned long long)v[i-1] >> (32-s));
  89    vn[0] = v[0] << s;
  90
  91    un = (unsigned *)alloca(4*(m + 1));
  92    un[m] = (unsigned long long)u[m-1] >> (32-s);
  93    for (i = m - 1; i > 0; i--)
  94       un[i] = (u[i] << s) | ((unsigned long long)u[i-1] >> (32-s));
  95    un[0] = u[0] << s;
  96
  97    for (j = m - n; j >= 0; j--) {       // Main loop.
  98       // Compute estimate qhat of q[j].
  99       qhat = (un[j+n]*b + un[j+n-1])/vn[n-1];
 100       rhat = (un[j+n]*b + un[j+n-1]) - qhat*vn[n-1];
 101 again:
 102       if (qhat >= b || qhat*vn[n-2] > b*rhat + un[j+n-2])
 103       { qhat = qhat - 1;
 104         rhat = rhat + vn[n-1];
 105         if (rhat < b) goto again;
 106       }
 107
 108       // Multiply and subtract.
 109       k = 0;
 110       for (i = 0; i < n; i++) {
 111          p = qhat*vn[i];
 112          t = un[i+j] - k - (p & 0xFFFFFFFFLL);
 113          un[i+j] = t;
 114          k = (p >> 32) - (t >> 32);
 115       }
 116       t = un[j+n] - k;
 117       un[j+n] = t;
 118
 119       q[j] = qhat;              // Store quotient digit.
 120       if (t < 0) {              // If we subtracted too
 121          q[j] = q[j] - 1;       // much, add back.
 122          k = 0;
 123          for (i = 0; i < n; i++) {
 124             t = (unsigned long long)un[i+j] + vn[i] + k;
 125             un[i+j] = t;
 126             k = t >> 32;
 127          }
 128          un[j+n] = un[j+n] + k;
 129       }
 130    } // End j.
 131    // If the caller wants the remainder, unnormalize
 132    // it and pass it back.
 133    if (r != NULL) {
 134       for (i = 0; i < n-1; i++)
 135          r[i] = (un[i] >> s) | ((unsigned long long)un[i+1] << (32-s));
 136       r[n-1] = un[n-1] >> s;
 137    }
 138    return 0;
 139 }
 140
 141 int errors;
 142
 143 void check(unsigned q[], unsigned r[],
 144            unsigned u[], unsigned v[],
 145            int m, int n,
 146            unsigned cq[], unsigned cr[]) {
 147    int i, szq;
 148
 149    szq = max(m - n + 1, 1);
 150    for (i = 0; i < szq; i++) {
 151       if (q[i] != cq[i]) {
 152          errors = errors + 1;
 153          dumpit("Error, dividend u =", m, u);
 154          dumpit("       divisor  v =", n, v);
 155          dumpit("For quotient,  got:", m-n+1, q);
 156          dumpit("        Should get:", m-n+1, cq);
 157          return;
 158       }
 159    }
 160    for (i = 0; i < n; i++) {
 161       if (r[i] != cr[i]) {
 162          errors = errors + 1;
 163          dumpit("Error, dividend u =", m, u);
 164          dumpit("       divisor  v =", n, v);
 165          dumpit("For remainder, got:", n, r);
 166          dumpit("        Should get:", n, cr);
 167          return;
 168       }
 169    }
 170    return;
 171 }
 172
 173 int main() {
 174    static unsigned test[] = {
 175    // m, n, u...,          v...,          cq...,  cr....
 176       1, 1, 3,             0,             1,      1,            // Error, divide by 0.
 177       1, 2, 7,             1,3,           0,      7,0,          // Error, n > m.
 178       2, 2, 0,0,           1,0,           0,      0,0,          // Error, incorrect remainder cr.
 179       1, 1, 3,             2,             1,      1,
 180       1, 1, 3,             3,             1,      0,
 181       1, 1, 3,             4,             0,      3,
 182       1, 1, 0,             0xffffffff,    0,      0,
 183       1, 1, 0xffffffff,    1,             0xffffffff, 0,
 184       1, 1, 0xffffffff,    0xffffffff,    1,      0,
 185       1, 1, 0xffffffff,    3,             0x55555555, 0,
 186       2, 1, 0xffffffff,0xffffffff, 1,     0xffffffff,0xffffffff, 0,
 187       2, 1, 0xffffffff,0xffffffff, 0xffffffff,        1,1,    0,
 188       2, 1, 0xffffffff,0xfffffffe, 0xffffffff,        0xffffffff,0, 0xfffffffe,
 189       2, 1, 0x00005678,0x00001234, 0x00009abc,        0x1e1dba76,0, 0x6bd0,
 190       2, 2, 0,0,           0,1,           0,      0,0,
 191       2, 2, 0,7,           0,3,           2,      0,1,
 192       2, 2, 5,7,           0,3,           2,      5,1,
 193       2, 2, 0,6,           0,2,           3,      0,0,
 194       1, 1, 0x80000000,  0x40000001, 0x00000001, 0x3fffffff,
 195       2, 1, 0x00000000,0x80000000, 0x40000001, 0xfffffff8,0x00000001, 0x00000008,
 196       2, 2, 0x00000000,0x80000000, 0x00000001,0x40000000, 0x00000001, 0xffffffff,0x3fffffff,
 197       2, 2, 0x0000789a,0x0000bcde, 0x0000789a,0x0000bcde,          1,          0,0,
 198       2, 2, 0x0000789b,0x0000bcde, 0x0000789a,0x0000bcde,          1,          1,0,
 199       2, 2, 0x00007899,0x0000bcde, 0x0000789a,0x0000bcde,          0, 0x00007899,0x0000bcde,
 200       2, 2, 0x0000ffff,0x0000ffff, 0x0000ffff,0x0000ffff,          1,          0,0,
 201       2, 2, 0x0000ffff,0x0000ffff, 0x00000000,0x00000001, 0x0000ffff, 0x0000ffff,0,
 202       3, 2, 0x000089ab,0x00004567,0x00000123, 0x00000000,0x00000001,   0x00004567,0x00000123, 0x000089ab,0,
 203       3, 2, 0x00000000,0x0000fffe,0x00008000, 0x0000ffff,0x00008000,   0xffffffff,0x00000000, 0x0000ffff,0x00007fff, // Shows that first qhat can = b + 1.
 204       3, 3, 0x00000003,0x00000000,0x80000000, 0x00000001,0x00000000,0x20000000,   0x00000003, 0,0,0x20000000, // Adding back step req'd.
 205       3, 3, 0x00000003,0x00000000,0x00008000, 0x00000001,0x00000000,0x00002000,   0x00000003, 0,0,0x00002000, // Adding back step req'd.
 206       4, 3, 0,0,0x00008000,0x00007fff, 1,0,0x00008000,   0xfffe0000,0, 0x00020000,0xffffffff,0x00007fff,  // Add back req'd.
 207       4, 3, 0,0x0000fffe,0,0x00008000, 0x0000ffff,0,0x00008000, 0xffffffff,0, 0x0000ffff,0xffffffff,0x00007fff,  // Shows that mult-sub quantity cannot be treated as signed.
 208       4, 3, 0,0xfffffffe,0,0x80000000, 0x0000ffff,0,0x80000000, 0x00000000,1, 0x00000000,0xfffeffff,0x00000000,  // Shows that mult-sub quantity cannot be treated as signed.
 209       4, 3, 0,0xfffffffe,0,0x80000000, 0xffffffff,0,0x80000000, 0xffffffff,0, 0xffffffff,0xffffffff,0x7fffffff,  // Shows that mult-sub quantity cannot be treated as signed.
 210    };
 211    int i, n, m, ncases, f;
 212    unsigned q[10], r[10];
 213    unsigned *u, *v, *cq, *cr;
 214
 215    printf("divmnu:\n");
 216    i = 0;
 217    ncases = 0;
 218    while (i < sizeof(test)/4) {
 219       m = test[i];
 220       n = test[i+1];
 221       u = &test[i+2];
 222       v = &test[i+2+m];
 223       cq = &test[i+2+m+n];
 224       cr = &test[i+2+m+n+max(m-n+1, 1)];
 225
 226       f = divmnu(q, r, u, v, m, n);
 227       if (f) {
 228          dumpit("Error return code for dividend u =", m, u);
 229          dumpit("                      divisor  v =", n, v);
 230          errors = errors + 1;
 231       }
 232       else
 233          check(q, r, u, v, m, n, cq, cr);
 234       i = i + 2 + m + n + max(m-n+1, 1) + n;
 235       ncases = ncases + 1;
 236    }
 237
 238    printf("%d errors out of %d cases; there should be 3.\n", errors, ncases);
 239    return 0;
 240 }