openpower/sv/bitmanip/mulmnu.c

   1 /* from hacker's delight, originsl by hannah suarez (hcs0) */
   2 // Computes the m+n-halfword product of n halfwords x m halfwords, unsigned.
   3 // Max line length is 57, to fit in hacker.book.
   4 #include <stdio.h>
   5 #include <stdlib.h>     //To define "exit", req'd by XLC.
   6
   7 // w[0], u[0], and v[0] contain the LEAST significant halfwords.
   8 // (The halfwords are in little-endian order).
   9 // This is Knuth's Algorithm M from [Knuth Vol. 2 Third edition (1998)]
  10 // section 4.3.1.  Picture is:
  11 //                   u[m-1] ... u[1] u[0]
  12 //                 x v[n-1] ... v[1] v[0]
  13 //                   --------------------
  14 //        w[m+n-1] ............ w[1] w[0]
  15
  16 void mulmnu(unsigned short w[], unsigned short u[],
  17    unsigned short v[], int m, int n) {
  18
  19    unsigned int k, t;
  20    int i, j;
  21
  22    for (i = 0; i < m; i++)
  23       w[i] = 0;
  24
  25    for (j = 0; j < n; j++) {
  26       k = 0;
  27       unsigned short phi[2000];
  28       unsigned short plo[2000];
  29       for (i = 0; i < m; i++) {
  30          unsigned product = u[i]*v[j] + w[i + j];
  31          phi[i] = product>>16;
  32          plo[i] = product;
  33       }
  34       for (i = 0; i < m; i++) {
  35          t = (phi[i]<<16) | plo[i]  + k;
  36          w[i + j] = t;          // (I.e., t & 0xFFFF).
  37          k = t >> 16;
  38       }
  39       w[j + m] = k;
  40    }
  41    return;
  42 }
  43
  44 int errors;
  45
  46 void check(unsigned short result[], unsigned short u[],
  47    unsigned short v[], int m, int n, unsigned short correct[]) {
  48    int i, j;
  49
  50    for (i = 0; i < m + n; i++) {
  51       if (correct[i] != result[i]) {
  52          errors = errors + 1;
  53          printf("Error, m = %d, n = %d, u = ", m, n);
  54          for (j = 0; j < m; j++) printf(" %04x", u[j]);
  55          printf(" v =");
  56          for (j = 0; j < n; j++) printf(" %04x", v[j]);
  57          printf("\nShould get:");
  58          for (j = 0; j < n+m; j++) printf(" %04x", correct[j]);
  59          printf("\n       Got:");
  60          for (j = 0; j < n+m; j++) printf(" %04x", result[j]);
  61          printf("\n");
  62          break;
  63       }
  64    }
  65 }
  66
  67 int main() {
  68    static unsigned short test[] = {
  69    //  m, n, u ..., v ..., result.
  70       1, 1, 7, 3,                  21,0,
  71       1, 1,      2, 0xFFFF,        0xFFFE,0x0001, // 2*FFFF = 0001_FFFE.
  72       1, 1, 0xFFFF, 0xFFFF,        1,0xFFFE,
  73       1, 2, 7, 5, 6,               35,42,0,
  74       1, 2, 65000, 63000, 64000,   0xBDC0,0x8414,0xF7F5,
  75       1, 3, 65535, 31000, 32000, 33000, 0x86E8,0xFC17,0xFC17,0x80E7,
  76       2, 3, 400, 300, 500, 100, 200, 0x0D40,0xE633,0xADB2,0xEA61,0,
  77       2, 3, 400, 65535, 500, 100, 65534, 0x0D40,0x9A4F,0xFE70,0x01F5,0xFFFD,
  78       4, 4, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535,
  79                 1,     0,     0,     0, 65534, 65535, 65535, 65535,
  80    };
  81    int i, n, m, ncases;
  82    unsigned short result[10];
  83    unsigned short *u, *v;
  84
  85    printf("mulmnu:\n");
  86    i = 0;
  87    ncases = 0;
  88    while (i < sizeof(test)/2) {
  89       m = test[i];
  90       n = test[i+1];
  91       u = &test[i+2];
  92       v = &test[i+2+m];
  93       mulmnu(result, u, v, m, n);
  94       check (result, u, v, m, n, &test[i+2+m+n]);
  95       mulmnu(result, v, u, n, m);       // Interchange operands.
  96       check (result, v, u, n, m, &test[i+2+m+n]);
  97       i = i + 2 + 2*(m + n);
  98       ncases = ncases + 1;
  99    }
 100
 101    if (errors == 0)
 102       printf("Passed all %d cases.\n", ncases);
 103 }