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


/* 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 */