1 b3f61791 2004-03-21 devnull #include "os.h"
2 b3f61791 2004-03-21 devnull #include <mp.h>
3 b3f61791 2004-03-21 devnull #include "dat.h"
6 cbeb0b26 2006-04-01 devnull /* from knuth's 1969 seminumberical algorithms, pp 233-235 and pp 258-260 */
8 cbeb0b26 2006-04-01 devnull /* mpvecmul is an assembly language routine that performs the inner */
9 cbeb0b26 2006-04-01 devnull /* loop. */
11 cbeb0b26 2006-04-01 devnull /* the karatsuba trade off is set empiricly by measuring the algs on */
12 cbeb0b26 2006-04-01 devnull /* a 400 MHz Pentium II. */
15 cbeb0b26 2006-04-01 devnull /* karatsuba like (see knuth pg 258) */
16 cbeb0b26 2006-04-01 devnull /* prereq: p is already zeroed */
17 b3f61791 2004-03-21 devnull static void
18 b3f61791 2004-03-21 devnull mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
20 b3f61791 2004-03-21 devnull mpdigit *t, *u0, *u1, *v0, *v1, *u0v0, *u1v1, *res, *diffprod;
21 b3f61791 2004-03-21 devnull int u0len, u1len, v0len, v1len, reslen;
22 b3f61791 2004-03-21 devnull int sign, n;
24 cbeb0b26 2006-04-01 devnull /* divide each piece in half */
25 b3f61791 2004-03-21 devnull n = alen/2;
26 b3f61791 2004-03-21 devnull if(alen&1)
28 b3f61791 2004-03-21 devnull u0len = n;
29 b3f61791 2004-03-21 devnull u1len = alen-n;
30 b3f61791 2004-03-21 devnull if(blen > n){
31 b3f61791 2004-03-21 devnull v0len = n;
32 b3f61791 2004-03-21 devnull v1len = blen-n;
34 b3f61791 2004-03-21 devnull v0len = blen;
35 b3f61791 2004-03-21 devnull v1len = 0;
38 b3f61791 2004-03-21 devnull u1 = a + u0len;
40 b3f61791 2004-03-21 devnull v1 = b + v0len;
42 cbeb0b26 2006-04-01 devnull /* room for the partial products */
43 b3f61791 2004-03-21 devnull t = mallocz(Dbytes*5*(2*n+1), 1);
44 b3f61791 2004-03-21 devnull if(t == nil)
45 b3f61791 2004-03-21 devnull sysfatal("mpkaratsuba: %r");
46 b3f61791 2004-03-21 devnull u0v0 = t;
47 b3f61791 2004-03-21 devnull u1v1 = t + (2*n+1);
48 b3f61791 2004-03-21 devnull diffprod = t + 2*(2*n+1);
49 b3f61791 2004-03-21 devnull res = t + 3*(2*n+1);
50 b3f61791 2004-03-21 devnull reslen = 4*n+1;
52 cbeb0b26 2006-04-01 devnull /* t[0] = (u1-u0) */
53 b3f61791 2004-03-21 devnull sign = 1;
54 b3f61791 2004-03-21 devnull if(mpveccmp(u1, u1len, u0, u0len) < 0){
55 b3f61791 2004-03-21 devnull sign = -1;
56 b3f61791 2004-03-21 devnull mpvecsub(u0, u0len, u1, u1len, u0v0);
58 b3f61791 2004-03-21 devnull mpvecsub(u1, u1len, u0, u1len, u0v0);
60 cbeb0b26 2006-04-01 devnull /* t[1] = (v0-v1) */
61 b3f61791 2004-03-21 devnull if(mpveccmp(v0, v0len, v1, v1len) < 0){
62 b3f61791 2004-03-21 devnull sign *= -1;
63 b3f61791 2004-03-21 devnull mpvecsub(v1, v1len, v0, v1len, u1v1);
65 b3f61791 2004-03-21 devnull mpvecsub(v0, v0len, v1, v1len, u1v1);
67 cbeb0b26 2006-04-01 devnull /* t[4:5] = (u1-u0)*(v0-v1) */
68 b3f61791 2004-03-21 devnull mpvecmul(u0v0, u0len, u1v1, v0len, diffprod);
70 cbeb0b26 2006-04-01 devnull /* t[0:1] = u1*v1 */
71 b3f61791 2004-03-21 devnull memset(t, 0, 2*(2*n+1)*Dbytes);
72 b3f61791 2004-03-21 devnull if(v1len > 0)
73 b3f61791 2004-03-21 devnull mpvecmul(u1, u1len, v1, v1len, u1v1);
75 cbeb0b26 2006-04-01 devnull /* t[2:3] = u0v0 */
76 b3f61791 2004-03-21 devnull mpvecmul(u0, u0len, v0, v0len, u0v0);
78 cbeb0b26 2006-04-01 devnull /* res = u0*v0<<n + u0*v0 */
79 b3f61791 2004-03-21 devnull mpvecadd(res, reslen, u0v0, u0len+v0len, res);
80 b3f61791 2004-03-21 devnull mpvecadd(res+n, reslen-n, u0v0, u0len+v0len, res+n);
82 cbeb0b26 2006-04-01 devnull /* res += u1*v1<<n + u1*v1<<2*n */
83 b3f61791 2004-03-21 devnull if(v1len > 0){
84 b3f61791 2004-03-21 devnull mpvecadd(res+n, reslen-n, u1v1, u1len+v1len, res+n);
85 b3f61791 2004-03-21 devnull mpvecadd(res+2*n, reslen-2*n, u1v1, u1len+v1len, res+2*n);
88 cbeb0b26 2006-04-01 devnull /* res += (u1-u0)*(v0-v1)<<n */
89 b3f61791 2004-03-21 devnull if(sign < 0)
90 b3f61791 2004-03-21 devnull mpvecsub(res+n, reslen-n, diffprod, u0len+v0len, res+n);
92 b3f61791 2004-03-21 devnull mpvecadd(res+n, reslen-n, diffprod, u0len+v0len, res+n);
93 b3f61791 2004-03-21 devnull memmove(p, res, (alen+blen)*Dbytes);
98 b3f61791 2004-03-21 devnull #define KARATSUBAMIN 32
101 b3f61791 2004-03-21 devnull mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
104 b3f61791 2004-03-21 devnull mpdigit d;
105 b3f61791 2004-03-21 devnull mpdigit *t;
107 cbeb0b26 2006-04-01 devnull /* both mpvecdigmuladd and karatsuba are fastest when a is the longer vector */
108 b3f61791 2004-03-21 devnull if(alen < blen){
109 b3f61791 2004-03-21 devnull i = alen;
110 b3f61791 2004-03-21 devnull alen = blen;
111 b3f61791 2004-03-21 devnull blen = i;
116 b3f61791 2004-03-21 devnull if(blen == 0){
117 b3f61791 2004-03-21 devnull memset(p, 0, Dbytes*(alen+blen));
121 b3f61791 2004-03-21 devnull if(alen >= KARATSUBAMIN && blen > 1){
122 cbeb0b26 2006-04-01 devnull /* O(n^1.585) */
123 b3f61791 2004-03-21 devnull mpkaratsuba(a, alen, b, blen, p);
124 b3f61791 2004-03-21 devnull } else {
125 cbeb0b26 2006-04-01 devnull /* O(n^2) */
126 b3f61791 2004-03-21 devnull for(i = 0; i < blen; i++){
127 b3f61791 2004-03-21 devnull d = b[i];
128 b3f61791 2004-03-21 devnull if(d != 0)
129 b3f61791 2004-03-21 devnull mpvecdigmuladd(a, alen, d, &p[i]);
135 b3f61791 2004-03-21 devnull mpmul(mpint *b1, mpint *b2, mpint *prod)
137 b3f61791 2004-03-21 devnull mpint *oprod;
139 b3f61791 2004-03-21 devnull oprod = nil;
140 b3f61791 2004-03-21 devnull if(prod == b1 || prod == b2){
141 b3f61791 2004-03-21 devnull oprod = prod;
142 b3f61791 2004-03-21 devnull prod = mpnew(0);
145 b3f61791 2004-03-21 devnull prod->top = 0;
146 b3f61791 2004-03-21 devnull mpbits(prod, (b1->top+b2->top+1)*Dbits);
147 b3f61791 2004-03-21 devnull mpvecmul(b1->p, b1->top, b2->p, b2->top, prod->p);
148 b3f61791 2004-03-21 devnull prod->top = b1->top+b2->top+1;
149 b3f61791 2004-03-21 devnull prod->sign = b1->sign*b2->sign;
150 b3f61791 2004-03-21 devnull mpnorm(prod);
152 b3f61791 2004-03-21 devnull if(oprod != nil){
153 b3f61791 2004-03-21 devnull mpassign(prod, oprod);
154 b3f61791 2004-03-21 devnull mpfree(prod);