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
mp_mult(MINT * a,MINT * b,MINT * c)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
m_mult(MINT * a,MINT * b,MINT * c)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 }
131