diff options
Diffstat (limited to 'src/libmp/port/mpmul.c')
| -rw-r--r-- | src/libmp/port/mpmul.c | 50 |
1 files changed, 25 insertions, 25 deletions
diff --git a/src/libmp/port/mpmul.c b/src/libmp/port/mpmul.c index dedd474a..a6ce9a0b 100644 --- a/src/libmp/port/mpmul.c +++ b/src/libmp/port/mpmul.c @@ -2,18 +2,18 @@ #include <mp.h> #include "dat.h" -// -// from knuth's 1969 seminumberical algorithms, pp 233-235 and pp 258-260 -// -// mpvecmul is an assembly language routine that performs the inner -// loop. -// -// the karatsuba trade off is set empiricly by measuring the algs on -// a 400 MHz Pentium II. -// - -// karatsuba like (see knuth pg 258) -// prereq: p is already zeroed +/* */ +/* from knuth's 1969 seminumberical algorithms, pp 233-235 and pp 258-260 */ +/* */ +/* mpvecmul is an assembly language routine that performs the inner */ +/* loop. */ +/* */ +/* the karatsuba trade off is set empiricly by measuring the algs on */ +/* a 400 MHz Pentium II. */ +/* */ + +/* karatsuba like (see knuth pg 258) */ +/* prereq: p is already zeroed */ static void mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p) { @@ -21,7 +21,7 @@ mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p) int u0len, u1len, v0len, v1len, reslen; int sign, n; - // divide each piece in half + /* divide each piece in half */ n = alen/2; if(alen&1) n++; @@ -39,7 +39,7 @@ mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p) v0 = b; v1 = b + v0len; - // room for the partial products + /* room for the partial products */ t = mallocz(Dbytes*5*(2*n+1), 1); if(t == nil) sysfatal("mpkaratsuba: %r"); @@ -49,7 +49,7 @@ mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p) res = t + 3*(2*n+1); reslen = 4*n+1; - // t[0] = (u1-u0) + /* t[0] = (u1-u0) */ sign = 1; if(mpveccmp(u1, u1len, u0, u0len) < 0){ sign = -1; @@ -57,35 +57,35 @@ mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p) } else mpvecsub(u1, u1len, u0, u1len, u0v0); - // t[1] = (v0-v1) + /* t[1] = (v0-v1) */ if(mpveccmp(v0, v0len, v1, v1len) < 0){ sign *= -1; mpvecsub(v1, v1len, v0, v1len, u1v1); } else mpvecsub(v0, v0len, v1, v1len, u1v1); - // t[4:5] = (u1-u0)*(v0-v1) + /* t[4:5] = (u1-u0)*(v0-v1) */ mpvecmul(u0v0, u0len, u1v1, v0len, diffprod); - // t[0:1] = u1*v1 + /* t[0:1] = u1*v1 */ memset(t, 0, 2*(2*n+1)*Dbytes); if(v1len > 0) mpvecmul(u1, u1len, v1, v1len, u1v1); - // t[2:3] = u0v0 + /* t[2:3] = u0v0 */ mpvecmul(u0, u0len, v0, v0len, u0v0); - // res = u0*v0<<n + u0*v0 + /* res = u0*v0<<n + u0*v0 */ mpvecadd(res, reslen, u0v0, u0len+v0len, res); mpvecadd(res+n, reslen-n, u0v0, u0len+v0len, res+n); - // res += u1*v1<<n + u1*v1<<2*n + /* res += u1*v1<<n + u1*v1<<2*n */ if(v1len > 0){ mpvecadd(res+n, reslen-n, u1v1, u1len+v1len, res+n); mpvecadd(res+2*n, reslen-2*n, u1v1, u1len+v1len, res+2*n); } - // res += (u1-u0)*(v0-v1)<<n + /* res += (u1-u0)*(v0-v1)<<n */ if(sign < 0) mpvecsub(res+n, reslen-n, diffprod, u0len+v0len, res+n); else @@ -104,7 +104,7 @@ mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p) mpdigit d; mpdigit *t; - // both mpvecdigmuladd and karatsuba are fastest when a is the longer vector + /* both mpvecdigmuladd and karatsuba are fastest when a is the longer vector */ if(alen < blen){ i = alen; alen = blen; @@ -119,10 +119,10 @@ mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p) } if(alen >= KARATSUBAMIN && blen > 1){ - // O(n^1.585) + /* O(n^1.585) */ mpkaratsuba(a, alen, b, blen, p); } else { - // O(n^2) + /* O(n^2) */ for(i = 0; i < blen; i++){ d = b[i]; if(d != 0) |