idigitkara.c

Arithmetic for integers of almost unlimited size for C and C++. Developed and copyrighted by Ra

/* Integer Version 2.0, RD, 20.1.93	idigitkara.c	*/
/*	use DigitVecMult, RD, 24.2.93			*/

/*	Karatsuba method for multiplication of long integers.	*/
/*	Parameters optimized for Sparc Version.		*/


#include <iint.h>
#include <idigit.h>

/*#ifndef KARATSUBA_LIMIT2
#define KARATSUBA_LIMIT2 15
#endif
*/
#define KARATSUBA_LIMIT2 16


/* Formulae:
    a == B^n a1 + a0;
    b == B^n b1 + b0;
    a b == (B^(2n) + B^n) a1 b1 + B^n (a1 - a0) (b0 - b1) + (B^n + 1) a0 b0;

	-----------------
V	|     a1 b1 	|
	-----------------
		-----------------
IV		|     a1 b1	|
		-----------------
		-----------------
III		||a1-a0||b0-b1|	|
		-----------------
		-----------------
II		|     a0 b0	|
		-----------------
			-----------------
I			|     a0 b0	|
			-----------------
	---------------------------------
	|	    prod		|
	---------------------------------

		|	|	|	|
		n3	n2	n	0

*/


static  DigitType vecadd4carry (/*sum, a, b, c, d, carry, k*/);
static  int vecadd3subcarry (/*accu, a, b, c, d, carry, k*/);
static  DigitType vecaddDigitType (/*accu, a, d, k*/);
static  DigitType vecaddscarry (/*accu, a, d, k*/);

void DigitVecKaratsubaM (prod, a, b, tmp, n2)
	DigitType *prod, *a, *b, *tmp;
	int n2;

/* 	In:	a[n2], b[n2];	n2 == 2^k * n0, with n0 <= KARATSUBA_LIMIT2
	Out:	prod[2*n2] = a[n2] * b[n2];
	Local:	tmp[2*n2] for intermediate results
	for n2>KARATSUBA_LIMIT2:	karatsuba-recursion with k-1
	else:		standard multiplication
*/
{	int sign, i, n, n3;
	DigitType carry;
	int subcarry;

	if (n2<=KARATSUBA_LIMIT2) {
	    /* use DigitVecMult and DigitVecMultAdd */
	    /* clearly n2 >= 1 */
	    prod[n2]=DigitVecMult(prod, a, b[0], n2);
            for ( i = 1; i < n2; i++) 
                prod[i+n2] = DigitVecMultAdd(&prod[i], a, b[i], n2);
	    return;
	}

	n=n2/2;
	n3=n2+n;

	DigitVecKaratsubaM(&prod[n2], a, b, prod, n);
	/* prod:	|    a0*b0	|		|	*/

	if (DigitVecGt(a, &a[n], n)) {
		sign=1;
		DigitVecSub(tmp, a, &a[n], n);	/* a0 - a1 */
	} else {
		sign=0;
		DigitVecSub(tmp, &a[n], a, n);	/* a1 - a0 */
	}
	/* tmp:		|	|	|	||a0-a1||	*/
	if (DigitVecGt(&b[n], b, n)) {
		sign^=1;
		DigitVecSub(&tmp[n], &b[n], b, n);	/* b1 - b0 */
	} else {
		DigitVecSub(&tmp[n], b, &b[n], n);	/* b0 - b1 */
	}
	/* tmp:		|	|	||b1-b0|||a0-a1||	*/


	DigitVecKaratsubaM(&tmp[n2], tmp, &tmp[n], prod, n);
	/* tmp:		||b1-b0|*|a0-a1|||b1-b0|||a0-a1||	*/

	DigitVecKaratsubaM(tmp, &a[n], &b[n], prod, n);
	/* tmp:		||b1-b0|*|a0-a1||    a1*b1	|	*/

	/*
			---------------------------------
	prod		|     a0 b0	|		|
			---------------------------------

			---------------------------------
	tmp		||a1-a0||b0-b1|	|     a1 b1	|
			---------------------------------
	*/

	if (sign==0) {		/* Einfach addieren */
	    for (i=0; i<n; i++)
		prod[i]=prod[n2+i];			/* I */
	    /* I + II + III + IV */
	    carry=vecadd4carry(&prod[n], 
			&prod[n3], &prod[n2], &tmp[n2], tmp, 0, n);
	    /* carry + II + III + IV + V */
	    carry=vecadd4carry(&prod[n2], 
			&prod[n3], &tmp[n3], &tmp[n], tmp, carry, n);
	    /* carry + V */
	    vecaddDigitType(&prod[n3], &tmp[n], carry, n);
	} else {		/* III subtrahieren */
	    for (i=0; i<n; i++)
		prod[i]=prod[n2+i];			/* I */
	    /* I + II + IV - III */
	    subcarry=vecadd3subcarry(&prod[n], 
			&prod[n3], &prod[n2], tmp, &tmp[n2], 0, n);
	    /* carry + II + IV + V - III */
	    subcarry=vecadd3subcarry(&prod[n2], 
			&prod[n3], &tmp[n], tmp, &tmp[n3], subcarry, n);
	    /* subcarry + V */
	    vecaddscarry (&prod[n3], &tmp[n], subcarry, n);
	}
}


static  DigitType vecadd4carry (accu, a, b, c, d, carry, k)
	DigitType *accu, *a, *b, *c, *d, carry;
	int k;
	/* accu[k]=a[k]+b[k]+c[k]+d[k]+carry; return carry; */
{	DigitType aa, bb, cc, dd, sum;
	for ( ; k>0; k--) {
		aa=*a++;
		bb=*b++;
		cc=*c++;
		dd=*d++;
		sum = aa+carry;
		carry = (sum < aa);
		sum += bb;
		carry += (sum < bb);
		sum += cc;
		carry += (sum < cc);
		sum += dd;
		carry += (sum < dd);
		*accu++= sum;
	}
	return carry;
}

static  int vecadd3subcarry (accu, a, b, c, d, carry, k)
	DigitType *accu, *a, *b, *c, *d;
	int carry, k;
	/* carry >= -1 */
	/* accu[k]=a[k]+b[k]+c[k]-d[k]+carry; return carry; */
{	DigitType aa, bb, cc, dd, sum, tmp, subc;
	if (carry < 0) {
		subc = 1;
		carry = 0;
	} else {
		subc=0;
	}
	for ( ; k>0; k--) {
		aa=*a++;
		bb=*b++;
		cc=*c++;
		dd=*d++;
		sum = aa+carry;
		carry = (sum < aa);
		sum += bb;
		carry += (sum < bb);
		sum += cc;
		carry += (sum < cc);
		tmp = sum - subc;
		subc = (tmp > sum);
		sum = tmp - dd;
		subc += (sum > tmp);
		*accu++= sum;
	}
	return carry-subc;
}

static  DigitType vecaddDigitType (accu, a, d, k)
	DigitType *accu, *a, d;
	int k;
	/* accu[k]=a[k]+d; return carry; */
{	DigitType sum, carry=d;
	for ( ; k>0; k-- ) {
		sum = *a++ + carry;
		carry = (sum < carry);
		*accu++ = sum;
	}
	return carry;
}
/* ACHTUNG:	k==0  ==>  carry==d */

static  DigitType vecaddscarry (accu, a, scarry, k)
	DigitType *accu, *a;
	int scarry, k;
	/* 2^BitsPerDigit > scarry >=-1 */
	/* accu[k]=a[k]+scarry; return carry; */
{	DigitType sum, carry, tmp;
	if (scarry >= 0) {
	    carry = scarry;
	    for ( ; k>0; k-- ) {
		sum = *a++ + carry;
		carry = (sum < carry);
		*accu++ = sum;
	    }
	    return carry;
	} else {
	    carry = 1;
	    for ( ; k>0; k-- ) {
		tmp = *a++;
		sum = tmp - carry;
		carry = (sum > tmp);
		*accu++ = sum;
	    }
	    return -carry;
	}
}
/* ACHTUNG:	k==0  ==>  carry==scarry */