igcd.c

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

/* Removed the "else" in bgcdkernel, TP, 27.01.95*/
/* Integer Version 2.1, RD, 7.5.93	igcd.c	*/
/* removed "else" in Ibgcd, RD, 7.5.93		*/
/* used Digitbgcd in Ibgcd, RD, 7.5.93		*/
/* "bgcdkernel", RD, 11.5.93			*/
/* possible bug in bgcdkernel fixed, RD, 22.7.93 */
/* another bug in bgcdkernel fixed (reported by Thomas P.), RD, 25.8.93 */

/* Integer Version 2.0, RD, 21.1.93	igcd.c	*/
/* Division with remainder in Ibgcd only for different lengths,
		RD, 12.2.93			*/
/* added two subtractions in Idgcd, Idxgcd before a division with
   remainder is done, added new function Idxgcd_left, VM, 09.12.97 */

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

void 
Idgcd(d, a, b)
    Integer *d;
    const Integer *a, *b;
/* d = gcd(a, b);   Euklidischer Algorithmus, Division mit Rest */
{
    Integer aa, bb, q, r;
    register pInteger pa, pb, pq, pr, copy;
    if (Ieq0(a))
    {
	IasI(d, b);
	d->sign = PLUS;
	return;
    }
    if (Ieq0(b))
    {
	IasI(d, a);
	d->sign = PLUS;
	return;
    }
    pa = &aa;
    pb = &bb;
    pq = &q;
    pr = &r;
    cIasI(pa, a);
    cIasI(pb, b);
    cI(pq);
    cI(pr);
    pa->sign = PLUS;
    pb->sign = PLUS;

    if (IltI(pa, pb))
      {
	copy = pa;
        pa = pb;
        pb = copy;
      }

    while (TRUE)
    {
        IasImiI(pr, pa, pb);
	if (IgeI(pr, pb))
	  {
	    ImiasI(pr, pb);
	    if (IgeI(pr, pb))
	      uIdiv(pq, pr, pr, pb);
	  }
	if (Ieq0(pr))
	    break;
	copy = pa;
	pa = pb;
	pb = pr;
	pr = copy;
    }
    IasI(d, pb);
    dI(pa);
    dI(pb);
    dI(pq);
    dI(pr);
}				/* Idgcd */

static DigitType 
Digitbgcd(a, b)
    DigitType a, b;
{
    int commontwo = 0;
    while (!(a & 1) && !(b & 1))
    {
	commontwo++;
	a >>= 1;
	b >>= 1;
    }
    while (!(a & 1))
	a >>= 1;
    while (!(b & 1))
	b >>= 1;
    while (TRUE)
	if (a >= b)
	{
	    a -= b;
	    if (!a)
		break;
	    while (!(a & 1))
		a >>= 1;
	}
	else
	{
	    b -= a;
	    while (!(b & 1))
		b >>= 1;
	}
    return b << commontwo;
}				/* Digitbgcd */

#ifndef USE_OLD_BGCD

static int 
bgcdkernel(d, a, al, b, bl)
    DigitType *d;
    DigitType *a;
    int al;
    DigitType *b;
    int bl;
/*	In:		a[0..al-1], b[0..bl-1]
	Compute:	d[0..dl-1] = gcd ( a[0..al-1], b[0..bl-1] );
	Out:		d[0..dl-1]; return dl;
	Where:		d[] has to be long enough to hold the result,
			a[0..al-1], b[0..bl-1] are destroyed,
			if dl>0 then d[dl-1] != 0,
			leading 0's are allowed,
			a == b == d are allowed.
	Note:		d[dl ...] are not changed.
			For preformance reasons the effective lengths
			of a and b should be similar.
	Uses:		DigitVecGt,
			DigitVecSri,
			DigitVecCsubto.
*/

{
    int i, digitshift = 0, bitshift = 0, j, k;
    DigitType a0, b0, s, t;

    /* get effective length of a and b */
    while (al > 0 && a[al - 1] == 0)
	al--;
    while (bl > 0 && b[bl - 1] == 0)
	bl--;

    /* return if a or b are 0 */
    if (al == 0)
    {
	for (i = 0; i < bl; i++)
	    d[i] = b[i];
	return bl;
    }
    if (bl == 0)
    {
	for (i = 0; i < al; i++)
	    d[i] = a[i];
	return al;
    }

    /* remove common powers of 2 */
    while ((*a & 1) == 0 && (*b & 1) == 0)
    {
	if (*a == 0 && *b == 0)
	{
	    a++;
	    b++;
	    al--;
	    bl--;
	    digitshift++;
	}
	else
	{
	    a0 = *a;
	    b0 = *b;
	    while ((a0 & 1) == 0 && (b0 & 1) == 0)
	    {
		a0 >>= 1;
		b0 >>= 1;
		bitshift++;
	    }
	    DigitVecSri(a, al, bitshift);
	    DigitVecSri(b, bl, bitshift);
	    if (a[al - 1] == 0)
		al--;
	    if (b[bl - 1] == 0)
		bl--;
	    break;
	}
    }

    /* make a and b odd */
    if ((*a & 1) == 0)
    {
	while (*a == 0)
	{
	    a++;
	    al--;
	}
	a0 = *a;
	i = 0;
	while ((a0 & 1) == 0)
	{
	    a0 >>= 1;
	    i++;
	}
	DigitVecSri(a, al, i);
	if (a[al - 1] == 0)
	    al--;
    }
    if ((*b & 1) == 0)
    {
	while (*b == 0)
	{
	    b++;
	    bl--;
	}
	b0 = *b;
	i = 0;
	while ((b0 & 1) == 0)
	{
	    b0 >>= 1;
	    i++;
	}
	DigitVecSri(b, bl, i);
	if (b[bl - 1] == 0)
	    bl--;
    }

    /* check for special case */
    if (al == 1 && bl == 1)
    {
	*b = Digitbgcd(*a, *b);
	bl = 1;
	goto bgcdfinalshift;
    }

    /* main loop */
    while (1)
    {
	/* here a, b are odd and al, bl are their lengths */
	if (bl > al || ((al == bl) && DigitVecGt(b, a, bl)))
	{
	    /* now b > a */
	    DigitVecCsubto(b, a, al);
	    while (b[bl - 1] == 0)
	    {
		bl--;
		if (bl == 1 && al == 1)
		{
		    *b = Digitbgcd(*a, *b);
		    goto bgcdfinalshift;
		}
	    }
	    if (*b == 0)
	    {
		while (*b == 0)
		{
		    b++;
		    bl--;
		    if (bl == 1 && al == 1)
		    {
			*b = Digitbgcd(*a, *b);
			goto bgcdfinalshift;
		    }
		}
		b0 = *b;
		i = 0;
	    }
	    else
	    {
		b0 = *b >> 1;
		i = 1;
	    }
	    while ((b0 & 1) == 0)
	    {
		b0 >>= 1;
		i++;
	    }
	    if (i)
		DigitVecSri(b, bl, i);
	    if (b[bl - 1] == 0)
	    {
		bl--;
		if (bl == 1 && al == 1)
		{
		    *b = Digitbgcd(*a, *b);
		    goto bgcdfinalshift;
		}
	    }
	}
	else
	{
	    /* now b <= a */
	    DigitVecCsubto(a, b, bl);
	    while (al > 0 && a[al - 1] == 0)
	    {
		al--;
		if (bl == 1 && al == 1)
		{
		    *b = Digitbgcd(*a, *b);
		    goto bgcdfinalshift;
		}
	    }

	    if (al == 0)
		goto bgcdfinalshift;

	    if (*a == 0)
	    {
		while (*a == 0)
		{
		    a++;
		    al--;
		    if (bl == 1 && al == 1)
		    {
			*b = Digitbgcd(*a, *b);
			goto bgcdfinalshift;
		    }
		}
		a0 = *a;
		i = 0;
	    }
	    else
	    {
		a0 = *a >> 1;
		i = 1;
	    }
	    while ((a0 & 1) == 0)
	    {
		a0 >>= 1;
		i++;
	    }
	    if (i)
		DigitVecSri(a, al, i);
	    if (a[al - 1] == 0)
	    {
		al--;
		if (bl == 1 && al == 1)
		{
		    *b = Digitbgcd(*a, *b);
		    goto bgcdfinalshift;
		}
	    }
	}
    }				/* main loop */

    /* now the result is b[0..bl-1] * powers of 2 */
bgcdfinalshift:
    for (i = 0; i < digitshift; i++)
	d[i] = 0;
    if (bitshift == 0)
    {
	for (i = 0; i < bl; i++)
	    d[digitshift + i] = b[i];
	return bl + digitshift;
    }
    /* now actual shift */
    j = bitshift;
    k = BitsPerDigit - j;
    s = 0;
    for (i = 0; i < bl; i++)
    {
	t = b[i];
	s |= (t << j);
	d[digitshift + i] = s;
	s = t >> k;
    }
    if (s)
    {
	d[digitshift + bl] = s;
	return digitshift + bl + 1;
    }
    else
	return digitshift + bl;
}				/* bgcdkernel */

void 
Ibgcd(d, a, b)
    Integer *d;
    const Integer *a, *b;
/*
 * d = gcd(a, b);   binaerer Algorithmus,
 * verwendet bgcdkernel
 * Achtung:
 * Sind die Zahlen a, b sehr unterschiedlich gross, so ist
 * Ibgcd relativ langsam. Deshalb ist es sinnvoll, vorher
 * eine Division mit Rest durchzuf"uhren.
 *
 */
{
    int l;
    Integer aa, bb;
    DigitType *v;

    if (!a->length)
    {
	IasI(d, b);
	d->sign = PLUS;
	return;
    }
    if (!b->length)
    {
	IasI(d, a);
	d->sign = PLUS;
	return;
    }
    cI(&aa);
    cI(&bb);
    /* Nun zuerst eine Division mit Rest */
    if (IgtI(a, b))
    {
	if (a->length >= (b->length << 1))
	{
	    IasIreI(&aa, a, b);
	    if (!aa.length)
	    {
		IasI(d, b);
		d->sign = PLUS;
		dI(&aa);
		dI(&bb);
		return;
	    }
	    IasI(&bb, b);
	}
	else
	{
	    IasI(&aa, b);
	    IasI(&bb, a);
	}
    }
    else
    {
	if (b->length >= (a->length << 1))
	{
	    IasIreI(&aa, b, a);
	    if (!aa.length)
	    {
		IasI(d, a);
		d->sign = PLUS;
		dI(&aa);
		dI(&bb);
		return;
	    }
	    IasI(&bb, a);
	}
	else
	{
	    IasI(&bb, b);
	    IasI(&aa, a);
	}
    }
    /* Nun ist bb >= aa */

    l = bgcdkernel(bb.vec, aa.vec, aa.length, bb.vec, bb.length);

    v = d->vec;
    d->vec = bb.vec;
    bb.vec = v;

    d->length = l;
    l = d->maxlength;
    d->maxlength = bb.maxlength;
    bb.maxlength = l;

    d->sign = PLUS;

    dI(&aa);
    dI(&bb);
}				/* Ibgcd */

#else

void 
Ibgcd(d, a, b)
    Integer *d;
    const Integer *a, *b;
/*
 * d = gcd(a, b);   binaerer Algorithmus,
 * direkt auf irv.c aufgesetzt.
 * Achtung:
 * Sind die Zahlen a, b sehr unterschiedlich gross, so ist
 * Ibgcd relativ langsam. Deshalb ist es sinnvoll, vorher
 * eine Division mit Rest durchzuf"uhren.
 *
 */
{
    register int l, al, bl;
    DigitType toshift = 0;
    Integer aa, bb;
    register DigitType *lp, *avec, *bvec;
    if (!a->length)
    {
	IasI(d, b);
	d->sign = PLUS;
	return;
    }
    if (!b->length)
    {
	IasI(d, a);
	d->sign = PLUS;
	return;
    }
    cI(&aa);
    cI(&bb);
    /* Nun zuerst eine Division mit Rest */
    if (IgtI(a, b))
    {
	if (a->length >= (b->length << 1))
	{
	    IasIreI(&aa, a, b);
	    if (!aa.length)
	    {
		IasI(d, b);
		d->sign = PLUS;
		dI(&aa);
		dI(&bb);
		return;
	    }
	}
	else
	{
	    IasI(&aa, a);
	}
	IasI(&bb, b);
    }
    else
    {
	if (b->length >= (a->length << 1))
	{