shithub: femtolisp

ref: ab876b803ac6efb812fddc8a154202b80632eb5e
dir: /3rd/mp/mpmul.c/

#include "platform.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
static void
mpkaratsuba(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
{
	mpdigit *t, *u0, *u1, *v0, *v1, *u0v0, *u1v1, *res, *diffprod;
	int u0len, u1len, v0len, v1len, reslen;
	int sign, n;

	// divide each piece in half
	n = alen/2;
	if(alen&1)
		n++;
	u0len = n;
	u1len = alen-n;
	if(blen > n){
		v0len = n;
		v1len = blen-n;
	} else {
		v0len = blen;
		v1len = 0;
	}
	u0 = a;
	u1 = a + u0len;
	v0 = b;
	v1 = b + v0len;

	// room for the partial products
	t = calloc(1, Dbytes*5*(2*n+1));
	if(t == nil)
		sysfatal("mpkaratsuba: %r");
	u0v0 = t;
	u1v1 = t + (2*n+1);
	diffprod = t + 2*(2*n+1);
	res = t + 3*(2*n+1);
	reslen = 4*n+1;

	// t[0] = (u1-u0)
	sign = 1;
	if(mpveccmp(u1, u1len, u0, u0len) < 0){
		sign = -1;
		mpvecsub(u0, u0len, u1, u1len, u0v0);
	} else
		mpvecsub(u1, u1len, u0, u1len, u0v0);

	// 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)
	mpvecmul(u0v0, u0len, u1v1, v0len, diffprod);

	// 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
	mpvecmul(u0, u0len, v0, v0len, u0v0);

	// 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
	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
	if(sign < 0)
		mpvecsub(res+n, reslen-n, diffprod, u0len+v0len, res+n);
	else
		mpvecadd(res+n, reslen-n, diffprod, u0len+v0len, res+n);
	memmove(p, res, (alen+blen)*Dbytes);

	free(t);
}

#define KARATSUBAMIN 32

void
mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
{
	int i;
	mpdigit d;
	mpdigit *t;

	// both mpvecdigmuladd and karatsuba are fastest when a is the longer vector
	if(alen < blen){
		i = alen;
		alen = blen;
		blen = i;
		t = a;
		a = b;
		b = t;
	}

	if(alen >= KARATSUBAMIN && blen > 1){
		// O(n^1.585)
		mpkaratsuba(a, alen, b, blen, p);
	} else {
		// O(n^2)
		for(i = 0; i < blen; i++){
			d = b[i];
			if(d != 0)
				mpvecdigmuladd(a, alen, d, &p[i]);
		}
	}
}

void
mpvectsmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p)
{
	int i;
	mpdigit *t;

	if(alen < blen){
		i = alen;
		alen = blen;
		blen = i;
		t = a;
		a = b;
		b = t;
	}
	if(blen == 0)
		return;
	for(i = 0; i < blen; i++)
		mpvecdigmuladd(a, alen, b[i], &p[i]);
}

void
mpmul(mpint *b1, mpint *b2, mpint *prod)
{
	mpint *oprod;

	oprod = prod;
	if(prod == b1 || prod == b2){
		prod = mpnew(0);
		prod->flags = oprod->flags;
	}
	prod->flags |= (b1->flags | b2->flags) & MPtimesafe;

	prod->top = 0;
	mpbits(prod, (b1->top+b2->top+1)*Dbits);
	if(prod->flags & MPtimesafe)
		mpvectsmul(b1->p, b1->top, b2->p, b2->top, prod->p);
	else
		mpvecmul(b1->p, b1->top, b2->p, b2->top, prod->p);
	prod->top = b1->top+b2->top+1;
	prod->sign = b1->sign*b2->sign;
	mpnorm(prod);

	if(oprod != prod){
		mpassign(prod, oprod);
		mpfree(prod);
	}
}