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