openpower/sv/biginteger/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[], unsigned short v[], int m,
  17             int n)
  18 {
  19
  20     unsigned int k, t;
  21     int i, j;
  22
  23     for (i = 0; i < m; i++)
  24         w[i] = 0;
  25
  26     for (j = 0; j < n; j++)
  27     {
  28         k = 0;
  29         unsigned short phi[2000];
  30         unsigned short plo[2000];
  31         for (i = 0; i < m; i++)
  32         {
  33             unsigned product = (unsigned)u[i] * v[j] + w[i + j];
  34             phi[i] = product >> 16;
  35             plo[i] = product;
  36         }
  37         for (i = 0; i < m; i++)
  38         {
  39             t = (((unsigned)phi[i] << 16) | plo[i]) + k;
  40             w[i + j] = t; // (I.e., t & 0xFFFF).
  41             k = t >> 16;
  42         }
  43         w[j + m] = k;
  44     }
  45     return;
  46 }
  47
  48 int errors;
  49
  50 void check(unsigned short result[], unsigned short u[], unsigned short v[],
  51            int m, int n, unsigned short correct[])
  52 {
  53     int i, j;
  54
  55     for (i = 0; i < m + n; i++)
  56     {
  57         if (correct[i] != result[i])
  58         {
  59             errors = errors + 1;
  60             printf("Error, m = %d, n = %d, u = ", m, n);
  61             for (j = 0; j < m; j++)
  62                 printf(" %04x", u[j]);
  63             printf(" v =");
  64             for (j = 0; j < n; j++)
  65                 printf(" %04x", v[j]);
  66             printf("\nShould get:");
  67             for (j = 0; j < n + m; j++)
  68                 printf(" %04x", correct[j]);
  69             printf("\n       Got:");
  70             for (j = 0; j < n + m; j++)
  71                 printf(" %04x", result[j]);
  72             printf("\n");
  73             break;
  74         }
  75     }
  76 }
  77
  78 int main()
  79 {
  80     static unsigned short test[] = {
  81         // clang-format off
  82    //  m, n, u ..., v ..., result.
  83       1, 1, 7, 3,                  21,0,
  84       1, 1,      2, 0xFFFF,        0xFFFE,0x0001, // 2*FFFF = 0001_FFFE.
  85       1, 1, 0xFFFF, 0xFFFF,        1,0xFFFE,
  86       1, 2, 7, 5, 6,               35,42,0,
  87       1, 2, 65000, 63000, 64000,   0xBDC0,0x8414,0xF7F5,
  88       1, 3, 65535, 31000, 32000, 33000, 0x86E8,0xFC17,0xFC17,0x80E7,
  89       2, 3, 400, 300, 500, 100, 200, 0x0D40,0xE633,0xADB2,0xEA61,0,
  90       2, 3, 400, 65535, 500, 100, 65534, 0x0D40,0x9A4F,0xFE70,0x01F5,0xFFFD,
  91       4, 4, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535,
  92                 1,     0,     0,     0, 65534, 65535, 65535, 65535,
  93       2, 2, 0xFF00,0xFF00, 0xFF00,0xFF00, 0,0xfe01,0xfc02,0xfe02
  94         // clang-format on
  95     };
  96     int i, n, m, ncases;
  97     unsigned short result[10];
  98     unsigned short *u, *v;
  99
 100     printf("mulmnu:\n");
 101     i = 0;
 102     ncases = 0;
 103     while (i < sizeof(test) / 2)
 104     {
 105         m = test[i];
 106         n = test[i + 1];
 107         u = &test[i + 2];
 108         v = &test[i + 2 + m];
 109         mulmnu(result, u, v, m, n);
 110         check(result, u, v, m, n, &test[i + 2 + m + n]);
 111         mulmnu(result, v, u, n, m); // Interchange operands.
 112         check(result, v, u, n, m, &test[i + 2 + m + n]);
 113         i = i + 2 + 2 * (m + n);
 114         ncases = ncases + 1;
 115     }
 116
 117     if (errors == 0)
 118         printf("Passed all %d cases.\n", ncases);
 119 }