6 /* from knuth's 1969 seminumberical algorithms, pp 233-235 and pp 258-260 */
8 /* mpvecmul is an assembly language routine that performs the inner */
11 /* the karatsuba trade off is set empiricly by measuring the algs on */
12 /* a 400 MHz Pentium II. */
15 /* karatsuba like (see knuth pg 258) */
16 /* prereq: p is already zeroed */
18 mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
20 mpdigit *t, *u0, *u1, *v0, *v1, *u0v0, *u1v1, *res, *diffprod;
21 int u0len, u1len, v0len, v1len, reslen;
24 /* divide each piece in half */
42 /* room for the partial products */
43 t = mallocz(Dbytes*5*(2*n+1), 1);
45 sysfatal("mpkaratsuba: %r");
48 diffprod = t + 2*(2*n+1);
54 if(mpveccmp(u1, u1len, u0, u0len) < 0){
56 mpvecsub(u0, u0len, u1, u1len, u0v0);
58 mpvecsub(u1, u1len, u0, u1len, u0v0);
61 if(mpveccmp(v0, v0len, v1, v1len) < 0){
63 mpvecsub(v1, v1len, v0, v1len, u1v1);
65 mpvecsub(v0, v0len, v1, v1len, u1v1);
67 /* t[4:5] = (u1-u0)*(v0-v1) */
68 mpvecmul(u0v0, u0len, u1v1, v0len, diffprod);
71 memset(t, 0, 2*(2*n+1)*Dbytes);
73 mpvecmul(u1, u1len, v1, v1len, u1v1);
76 mpvecmul(u0, u0len, v0, v0len, u0v0);
78 /* res = u0*v0<<n + u0*v0 */
79 mpvecadd(res, reslen, u0v0, u0len+v0len, res);
80 mpvecadd(res+n, reslen-n, u0v0, u0len+v0len, res+n);
82 /* res += u1*v1<<n + u1*v1<<2*n */
84 mpvecadd(res+n, reslen-n, u1v1, u1len+v1len, res+n);
85 mpvecadd(res+2*n, reslen-2*n, u1v1, u1len+v1len, res+2*n);
88 /* res += (u1-u0)*(v0-v1)<<n */
90 mpvecsub(res+n, reslen-n, diffprod, u0len+v0len, res+n);
92 mpvecadd(res+n, reslen-n, diffprod, u0len+v0len, res+n);
93 memmove(p, res, (alen+blen)*Dbytes);
98 #define KARATSUBAMIN 32
101 mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
107 /* both mpvecdigmuladd and karatsuba are fastest when a is the longer vector */
117 memset(p, 0, Dbytes*(alen+blen));
121 if(alen >= KARATSUBAMIN && blen > 1){
123 mpkaratsuba(a, alen, b, blen, p);
126 for(i = 0; i < blen; i++){
129 mpvecdigmuladd(a, alen, d, &p[i]);
135 mpmul(mpint *b1, mpint *b2, mpint *prod)
140 if(prod == b1 || prod == b2){
146 mpbits(prod, (b1->top+b2->top+1)*Dbits);
147 mpvecmul(b1->p, b1->top, b2->p, b2->top, prod->p);
148 prod->top = b1->top+b2->top+1;
149 prod->sign = b1->sign*b2->sign;
153 mpassign(prod, oprod);