idigitvec.c

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

/* Corrected DigitVecSri, TP, 27.01.95 */
/* DigitVecAddCarry and DigitVecSubCarry removed, RD, 14.7.93 */
/* Integer Version 2.1, DigitVecMultSub changed, RD, 30.4.93	*/
/* Integer Version 2.0, RD, 18.1.93	idigitvec.c	*/
/*	Sparc Version 8, RD, 11.2.93	*/
/* DigitVecSri: try to avoid memory access delay, RD, 15.2.93 */

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

/*******************************
	First the standard DigitVec operations 
*********************************/

DigitType  DigitVecAdd (sum, a, b, l)
	DigitType *sum, *a, *b;
	int l;
     /*	sum[0..l-1] = a[0..l-1] + b[0..l-1]; return CARRY;   */
     /* This notation means:
	Add the integers represented by the vectors a and b,
	write the lowest l digits of the result into the vector sum
	and return the carry of this operation.
     */
{	DigitType accu, carry, tmp;
	carry=0;
	for ( ; l>0; l-- ) {
		accu = *a++ + carry;
		carry = (accu < carry);
		tmp = *b++;
		accu += tmp;
		carry += (accu < tmp);
		*sum++ = accu;
	}
	return carry;
}		/* DigitVecAdd */

DigitType  DigitVecSub (diff, a, b, l)
	DigitType *diff, *a, *b;
	int l;
     /*	diff[0..l-1] = a[0..l-1] - b[0..l-1]; return CARRY;   */
{	DigitType accu, carry, tmp;
	carry=0;
	for ( ; l>0; l-- ) {
		tmp = *a++;
		accu = tmp - carry;
		carry = (accu > tmp);
		tmp = accu - *b++;
		carry += (tmp > accu);
		*diff++ = tmp;
	}
	return carry;
}		/* DigitVecSub */

#define DDH (BitsPerDigit/2)

#ifdef USE_C_DIGITVECMULT

DigitType  DigitVecMult (res, a, m, l)
	DigitType *res, *a, m;
	int l;
     /*	res[0..l-1] = a[0..l-1]*m; return CARRY;   */
{       DigitType carry;
        DigitType aa, low1, low2, tmp1, tmp2, tmp3;
        DigitType ml, mh, al, ah;

        ml = (m << DDH) >> DDH;
        mh = m >> DDH;
	carry=0;

        for ( ; l>0; l--) {
            aa = *a++;
            al = (aa << DDH) >> DDH;
            ah = aa >> DDH;

            tmp1 = ml * al;
            low2 = carry + tmp1;
            carry = (low2 < carry);

            tmp2 = ml * ah;
            low1 = low2 + (tmp2 << DDH);
            carry += (low1 < low2) + (tmp2 >> DDH);

            tmp3 = mh * al;
            low2 = low1 + (tmp3 << DDH);
            carry += (low2 < low1) + (tmp3 >> DDH);

            *res++ = low2;
            carry += mh * ah;
        }
        return carry;
}               /* DigitVecMult */

#endif

#ifdef USE_C_DIGITVECMULTADD

DigitType  DigitVecMultAdd (res, a, m, l)
	DigitType *res, *a, m;
	int l;
     /*	res[0..l-1] += a[0..l-1]*m; return CARRY;   */
     /* Special function for multiple digit multiplication */
{       DigitType carry;
        DigitType aa, low1, low2, tmp1, tmp2, tmp3;
        DigitType ml, mh, al, ah;

        ml = (m << DDH) >> DDH;
        mh = m >> DDH;
	carry=0;

        for ( ; l>0; l--) {
            aa = *a++;
            al = (aa << DDH) >> DDH;
            ah = aa >> DDH;

            low1  = carry + *res;
            carry = (low1 < carry);

            tmp1 = ml * al;
            low2 = low1 + tmp1;
            carry += (low2 < low1);

            tmp2 = ml * ah;
            low1 = low2 + (tmp2 << DDH);
            carry += (low1 < low2) + (tmp2 >> DDH);

            tmp3 = mh * al;
            low2 = low1 + (tmp3 << DDH);
            carry += (low2 < low1) + (tmp3 >> DDH);

            *res++ = low2;
            carry += mh * ah;
        }
        return carry;
}               /* DigitVecMultAdd */

#endif

#ifdef USE_C_DIGITVECMULTSUB

DigitType  DigitVecMultSub (res, a, m, l)
	DigitType *res, *a, m;
	int l;
     /*	res[0..l-1] -= a[0..l-1]*m; return CARRY;   */
     /* Special function for multiple digit division */
{       DigitType carry;
        DigitType aa, low1, low2, tmp1, tmp2, tmp3;
        DigitType ml, mh, al, ah;

        ml = (m << DDH) >> DDH;
        mh = m >> DDH;
	carry=0;

        for ( ; l>0; l--) {
            aa = *a++;
            al = (aa << DDH) >> DDH;
            ah = aa >> DDH;

            tmp1 = ml * al;
            low2 = carry + tmp1;
            carry = (low2 < carry);

            tmp2 = ml * ah;
            low1 = low2 + (tmp2 << DDH);
            carry += (low1 < low2) + (tmp2 >> DDH);

            tmp3 = mh * al;
            low2 = low1 + (tmp3 << DDH);
            carry += (low2 < low1) + (tmp3 >> DDH);

	    tmp1 = *res;
	    low1 = tmp1 - low2;
	    carry += (low1 > tmp1);

            *res++ = low1;
            carry += mh * ah;
        }
        return carry;
}               /* DigitVecMultSub */

#endif

DigitType  DigitVecDiv (quot, a, d, l)
	DigitType *quot, *a, d;
	int l;
     /*	quot[0..l-1] = a[0..l-1]/m; 
	return a[0..l-1]%m;
     */
{       DigitType carry=0;
        for (quot+=l, a+=l; l>0; l--)
		carry=DigitDiv(--quot, carry, *--a, d);
        return carry;
}		/* DigitVec_div */


/*******************************
	Now some special DigitVec operations 
*********************************/

void DigitVecCsubto (a, b, l)
	DigitType *a, *b;
	int l;
        /* a[]-=b[l]; */
{	DigitType carry, x, y;
	carry=0;
        for ( ; l>0; l--) {
		x = *a;
		y = x - carry;
		carry = (y > x);
		x = y - *b++;
		carry += (x > y);
		*a++ = x;
	}
        while (carry) {
		x = *a;
		y = x - carry;
		carry = (y > x);
		*a++ = y;
	}
}		/* DigitVecCsubto */

int DigitVecCadd (sum, a, b, la, lb)
	DigitType *sum, *a, *b;
	int la, lb;
	/* sum[]=a[la]+b[lb]; return sum->length */
{       DigitType carry=0, accu, tmp;
	int countmax;
	if (la>lb) {
		countmax=la;
		la-=lb;
	        for ( ; lb>0; lb--) {
			accu = *a++ + carry;
			carry = (accu < carry);
			tmp = *b++;
			accu += tmp;
			carry += (accu < tmp);
			*sum++ = accu;
		}
	        for ( ; la>0; la--) {
			accu = *a++ + carry;
			carry = (accu < carry);
			*sum++ = accu;
		}
	} else {
		countmax=lb;
		lb-=la;
	        for ( ; la>0; la--) {
			accu = *a++ + carry;
			carry = (accu < carry);
			tmp = *b++;
			accu += tmp;
			carry += (accu < tmp);
			*sum++ = accu;
		}
	        for ( ; lb>0; lb--) {
			accu = *b++ + carry;
			carry = (accu < carry);
			*sum++ = accu;
    	}	}
	*sum = carry;
	if (carry)
		return countmax+1;
	else
		return countmax;
}	/* DigitVecCadd */

int DigitVecCsub (diff, a, b, la, lb)
	DigitType *diff, *a, *b;
	int la, lb;
	/* diff[]=a[la]-b[lb]; return diff->length */
{       DigitType carry=0, x, y;
	int l=la-lb;
	for ( ; lb>0; lb--) {
		x = *a++;
		y = x - carry;
		carry = (y > x);
		x = y - *b++;
		carry += (x > y);
		*diff++ = x;
	}
	while (carry) {
		x = *a++;
		y = x - carry;
		carry = (y > x);
		*diff++ = y;
		l--;
	}
	for ( ; l>0; l--) {
		*diff++=*a++;
	}
	diff--;
	while ((la>0)&&(! *diff)) {
		diff--;
		la--;
	}
	return la;
}	/* DigitVecCsub */

BOOLEAN DigitVecSr1 (u, l) 
	DigitType *u;
	int l;
	/* b=u[l]%2; u[l]/=2; return b; */
{	DigitType b, bold, c;
	u=u + l;
	b=bold=0;
	while(l) {
		b=*--u;
		c=b>>1;
		c|=(bold<<(BitsPerDigit-1));
		bold=b;
		*u=c;
		l--;
	}
	b&=1;
	return (BOOLEAN) b;
}		/* DigitVecSr1 */

#ifdef USE_OLD_VECSRI

void DigitVecSri (u, l, i)
	DigitType *u;
	int l, i;
	/* b=u[l]%2^i; u[l]/=2^i; 	0<i<BitsPerDigit  */
{	DigitType b, bold, c;
	u=u + l;
	bold=0;
	while(l) {
		b=*--u;
		c=b>>i;
		c|=(bold<<(BitsPerDigit-i));
		bold=b;
		*u=c;
		l--;
	}
}		/* DigitVecSri */

#else

void DigitVecSri (u, l, i)
	DigitType *u;
	int l, i;
	/* b=u[l]%2^i; u[l]/=2^i; 	0<i<BitsPerDigit  */
{	DigitType b, bold, c;
	if ((!l) || (!i))
		return;
	u = u + l - 1;
	b = *u;
	bold=0;
	l--;
	while(l) {
		c=b>>i;
		c|=(bold<<(BitsPerDigit-i));
		*u=c;
		bold=b;
		b=*--u;
		l--;
	}
	c=b>>i;
	c|=(bold<<(BitsPerDigit-i));
	bold=b;
	*u=c;
}		/* DigitVecSri */

#endif


/*******************************
	Now some general non-critical DigitVec operations 
*********************************/

BOOLEAN DigitVecEq (a, b, l)
	DigitType *a, *b;
	int l;
        /* return a[l]==b[l]; */
{       for ( ; l>0; l--)
                if (*a++ != *b++)
                        return FALSE;
        return TRUE;
}

BOOLEAN DigitVecGt (a, b, l)
	DigitType *a, *b;
	int l;
        /* return a[l]>b[l] lexikographisch */
{       for (a+=l, b+=l; l>0; l--) {
		DigitType aa, bb;
		aa=*--a;
		bb=*--b;
                if (aa > bb)
                        return TRUE;
                else if (aa < bb)
                        return FALSE;
	}
        return FALSE;
}