usr/src/lib/libmp/common/mult.c

   1 /*      Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
   2 /*        All Rights Reserved   */
   3
   4
   5 /*
   6  * Copyright (c) 1980 Regents of the University of California.
   7  * All rights reserved.  The Berkeley software License Agreement
   8  * specifies the terms and conditions for redistribution.
   9  */
  10 /*      Portions Copyright(c) 1988, Sun Microsystems Inc.       */
  11 /*      All Rights Reserved                                     */
  12
  13 /*
  14  * Copyright (c) 1997, by Sun Microsystems, Inc.
  15  * All rights reserved.
  16  */
  17
  18 #ident  "%Z%%M% %I%     %E% SMI"        /* SVr4.0 1.1   */
  19
  20 /* LINTLIBRARY */
  21
  22 #include <mp.h>
  23 #include <sys/types.h>
  24 #include "libmp.h"
  25
  26 static void m_mult(MINT *, MINT *, MINT *);
  27
  28 void
  29 mp_mult(MINT *a, MINT *b, MINT *c)
  30 {
  31         struct mint x, y;
  32         int sign;
  33
  34         _mp_mcan(a);
  35         _mp_mcan(b);
  36         if (a->len == 0 || b->len == 0) {
  37                 _mp_xfree(c);
  38                 return;
  39         }
  40         sign = 1;
  41         x.len = y.len = 0;
  42         _mp_move(a, &x);
  43         _mp_move(b, &y);
  44         if (a->len < 0) {
  45                 x.len = -x.len;
  46                 sign = -sign;
  47         }
  48         if (b->len < 0) {
  49                 y.len = -y.len;
  50                 sign = -sign;
  51         }
  52         _mp_xfree(c);
  53         if (x.len < y.len) {
  54                 m_mult(&x, &y, c);
  55         } else {
  56                 m_mult(&y, &x, c);
  57         }
  58         if (sign < 0)
  59                 c->len = -c->len;
  60         if (c->len == 0)
  61                 _mp_xfree(c);
  62         _mp_xfree(&x);
  63         _mp_xfree(&y);
  64 }
  65
  66 /*
  67  * Knuth  4.3.1, Algorithm M
  68  */
  69 static void
  70 m_mult(MINT *a, MINT *b, MINT *c)
  71 {
  72         int i, j;
  73         int sum;
  74         short bcache;
  75         short *aptr;
  76         short *bptr;
  77         short *cptr;
  78         short fifteen = 15;
  79         int alen;
  80         int blen;
  81
  82 #define BASEBITS        (8 * (unsigned int)sizeof (short) - 1)
  83 #define BASE            (1 << BASEBITS)
  84 #define LOWBITS         (BASE - 1)
  85
  86         alen = a->len;
  87         blen = b->len;
  88
  89         c->len = alen + blen;
  90         c->val = _mp_xalloc(c->len, "m_mult");
  91
  92         aptr = a->val;
  93         bptr = b->val;
  94         cptr = c->val;
  95
  96         sum = 0;
  97         bcache = *bptr++;
  98         for (i = alen; i > 0; i--) {
  99                 sum += *aptr++ * bcache;
 100                 *cptr++ = (short)(sum & LOWBITS);
 101                 if (sum >= 0)
 102                         sum >>= fifteen;
 103                 else
 104                         sum = 0xfffe0000 | (sum >> fifteen);
 105         }
 106         *cptr = (short)sum;
 107         aptr -= alen;
 108         cptr -= alen;
 109         cptr++;
 110
 111         for (j = blen - 1; j > 0; j--) {
 112                 sum = 0;
 113                 bcache = *bptr++;
 114                 for (i = alen; i > 0; i--) {
 115                         sum += *aptr++ * bcache + *cptr;
 116                         *cptr++ = (short)(sum & LOWBITS);
 117                         if (sum >= 0)
 118                                 sum >>= fifteen;
 119                         else
 120                                 sum = 0xfffe0000 | (sum >> fifteen);
 121                 }
 122                 *cptr = (short)sum;
 123                 aptr -= alen;
 124                 cptr -= alen;
 125                 cptr++;
 126         }
 127         if (c->val[c->len-1] == 0) {
 128                 c->len--;
 129         }
 130 }