igcd.c

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

	    IasIreI(&bb, b, a);
	    if (!bb.length)
	    {
		IasI(d, a);
		d->sign = PLUS;
		dI(&aa);
		dI(&bb);
		return;
	    }
	}
	else
	{
	    IasI(&bb, b);
	}
	IasI(&aa, a);
    }


    avec = aa.vec;
    bvec = bb.vec;
    al = aa.length;
    bl = bb.length;
    while (!(*avec & 1) && !(*bvec & 1))
    {
	toshift++;
	DigitVecSr1(avec, al);
	if (!avec[al - 1])
	    al--;
	DigitVecSr1(bvec, bl);
	if (!bvec[bl - 1])
	    bl--;
    }
    while (!(*avec & 1))
    {
	DigitVecSr1(avec, al);
	if (!avec[al - 1])
	    al--;
    }
    while (!(*bvec & 1))
    {
	DigitVecSr1(bvec, bl);
	if (!bvec[bl - 1])
	    bl--;
    }
    while (TRUE)
    {
	if (bl > al || ((al == bl) && DigitVecGt(bvec, avec, bl)))
	{
	    lp = avec;
	    avec = bvec;
	    bvec = lp;
	    l = al;
	    al = bl;
	    bl = l;
	}
	/* now  a >= b */
	lp = &(avec[al - 1]);
	DigitVecCsubto(avec, bvec, bl);
	if (!*lp)
	{
	    while ((al > 0) && (!*lp))
	    {
		al--;
		lp--;
	    }
	    if (!al)
		break;
	    if ((al == 1) && (bl == 1))
	    {
		*bvec = Digitbgcd(*bvec, *avec);
		break;
	    }
	}
	if (al)
	    while (!(*avec & 1))
	    {
		register DigitType low;
		l = 1;
		low = *avec >> 1;
		while (!(low & 1) && l < BitsPerDigit - 1)
		{
		    l++;
		    low >>= 1;
		}
		DigitVecSri(avec, al, l);
		if (!avec[al - 1])
		{
		    al--;
		    if ((al == 1) && (bl == 1))
		    {
			*bvec = Digitbgcd(*bvec, *avec);
			goto end;
		    }
		}
	    }
	else
	    break;		/* gcd in  bvec[0..bl-1] */
    }
end:if (avec == aa.vec)
    {
	bb.length = bl;
	bb.sign = PLUS;
	IasIslD(d, &bb, toshift);
	dI(&aa);
	dI(&bb);
	return;
    }
    else
    {
	aa.length = bl;
	aa.sign = PLUS;
	IasIslD(d, &aa, toshift);
	dI(&aa);
	dI(&bb);
	return;
    }
}				/* Ibgcd */

#endif

/*************************************************/

void 
Idxgcd(d, u, v, a, b)
    Integer *d, *u, *v;
    const Integer *a, *b;
/* d = gcd(a, b) = ua + vb;   Euklid-Lenstra-Berlekamp */
{
    Integer aa, bb, wua, wva, wub, wvb, q, r;
    register pInteger pa, pb, pq, pr, ua, va, ub, vb, copy;
    int exchanged;

    if (Ieq0(a))
    {
        IasI(d, b);
	Ias0(u);
	Ias1(v);
	v->sign = b->sign;
	d->sign = PLUS;
	return;
    }
    if (Ieq0(b))
    {
	IasI(d, a);
	Ias1(u);
	Ias0(v);
	u->sign = a->sign;
	d->sign = PLUS;
	return;
    }
    pa = &aa;
    pb = &bb;
    pq = &q;
    pr = &r;
    ua = &wua;
    va = &wva;
    ub = &wub;
    vb = &wvb;
    cIasI(pa, a);
    cIasI(pb, b);
    cI(pq);
    cI(pr);
    pa->sign = PLUS;
    pb->sign = PLUS;

    if (IgeI(pa, pb))
      {
	cIasint(ua, 1);
	ua->sign = a->sign;
	cIasint(va, 0);		/* pa == ua*a + va*b */
	cIasint(ub, 0);
	cIasint(vb, 1);
	vb->sign = b->sign;	/* pb == ub*a + vb*b */
        exchanged = 0;
      }
    else
      {
	copy = pa;
        pa = pb;
        pb = copy;

	cIasint(ua, 1);
	ua->sign = b->sign;
	cIasint(va, 0);		/* pa == ua*b + va*a */
	cIasint(ub, 0);
	cIasint(vb, 1);
	vb->sign = a->sign;	/* pb == ub*b + vb*a */
        exchanged = 1;
      }

    while (TRUE)
    {
       IasImiI(pr, pa, pb);
       if (IltI(pr, pb))
	 {
	     if (Ieq0(pr))
	       break;
	     copy = pa;
	     pa = pb;
	     pb = pr;
	     pr = copy;
	     copy = ua;
	     ua = ub;
	     ub = copy;
	     ImiasI(ub, ua);
	     copy = va;
	     va = vb;
	     vb = copy;
	     ImiasI(vb, va);
	 }
       else
	 { 
	   ImiasI(pr, pb);
	   if (Ieq0(pr))
	     break;
	   if (IltI(pr, pb))
	     Iasint(pq, 2);
	   else
	     {
	       uIdiv(pq, pr, pr, pb);
	       Iinc(pq);
	       Iinc(pq);
	     }
	   if (Ieq0(pr))
	     break;
	   copy = pa;
	   pa = pb;
	   pb = pr;
	   pr = copy;
	   IasImuI(pr, pq, ub);
	   copy = ua;
	   ua = ub;
	   ub = copy;
	   ImiasI(ub, pr);
	   IasImuI(pr, pq, vb);
	   copy = va;
	   va = vb;
	   vb = copy;
	   ImiasI(vb, pr);
	 }
    }
    IasI(d, pb);

    if (exchanged == 0)
      {
	IasI(u, ub);
	IasI(v, vb);
      }
    else
      {
	IasI(v, ub);
	IasI(u, vb);
      }

    dI(pa);
    dI(pb);
    dI(pq);
    dI(pr);
    dI(ua);
    dI(ub);
    dI(va);
    dI(vb);
}				/* Idxgcd */

void 
Idxgcd_left(d, u, a, b)
    Integer *d, *u;
    const Integer *a, *b;
/* d = gcd(a, b) = ua + ?*b;   Euklid-Lenstra-Berlekamp */
{
    Integer aa, bb, wua, wub, q, r;
    register pInteger pa, pb, pq, pr, ua, ub, copy;
    if (Ieq0(a))
    {
        IasI(d, b);
        Ias0(u);
        d->sign = PLUS;
        return;
    }
    if (Ieq0(b))
    {
        IasI(d, a);
        Ias1(u);
        u->sign = a->sign;
        d->sign = PLUS;
        return;
    }

    pa = &aa;
    pb = &bb;
    pq = &q;
    pr = &r;
    ua = &wua;
    ub = &wub;

    cIasI(pa, a);
    cIasI(pb, b);
    cI(pq);
    cI(pr);
    pa->sign = PLUS;
    pb->sign = PLUS;


    if (!IltI(pa, pb))
      {
	cIasint(ua, 1);
	ua->sign = a->sign;
	cIasint(ub, 0);
      }
    else
      {
	copy = pa;
        pa = pb;
        pb = copy;
	cIasint(ua, 0);
	cIasint(ub, 1);
	ub->sign = a->sign;
      }
    

    while (TRUE)
    {
      IasImiI(pr, pa, pb);
      if (IltI(pr, pb))
	{
	  if (Ieq0(pr))
	    break;
	  copy = pa;
	  pa = pb;
	  pb = pr;
	  pr = copy;
	  copy = ua;
	  ua = ub;
	  ub = copy;
	  ImiasI(ub, ua);
	}
      else
	{
	  ImiasI(pr, pb);
	  if (Ieq0(pr))
	    break;

	  if (IltI(pr, pb))
	    Iasint(pq, 2);
	  else
	    {
	      uIdiv(pq, pr, pr, pb);
	      Iinc(pq);
	      Iinc(pq);
	    }
	  if (Ieq0(pr))
	    break;
	  copy = pa;
	  pa = pb;
	  pb = pr;
	  pr = copy;
	  IasImuI(pr, pq, ub);
	  copy = ua;
	  ua = ub;
	  ub = copy;
	  ImiasI(ub, pr);
	}
    }
    IasI(d, pb);
    IasI(u, ub);
    dI(pa);
    dI(pb);
    dI(pq);
    dI(pr);
    dI(ua);
    dI(ub);
}              /* Idxgcd_left */
              

void 
Ibxgcd(d, u, v, a, b)
    Integer *d, *u, *v;
    const Integer *a, *b;
/* d = gcd(a, b) = u*a + v*b;   binaerer Algorithmus */
{
    DigitType toshift = 0;
    Integer aa, bb, A, B, wxa, wxb, wya, wyb;
    register pInteger pa, pb, pA, pB, xa, xb, ya, yb, swap;
    if (Ieq0(a))
    {
	IasI(d, b);
	Ias0(u);
	Ias1(v);
	v->sign = b->sign;
	d->sign = PLUS;
	return;
    }
    if (Ieq0(b))
    {
	IasI(d, a);
	Ias1(u);
	Ias0(v);
	u->sign = a->sign;
	d->sign = PLUS;
	return;
    }
    pa = &aa;
    pb = &bb;
    xa = &wxa;
    xb = &wxb;
    ya = &wya;
    yb = &wyb;
    pA = &A;
    pB = &B;
    cIasI(pa, a);
    cIasI(pb, b);
    while (Ieven(pa) && Ieven(pb))
    {
	toshift++;
	Isr1(pa);
	Isr1(pb);
    }
    cIasI(pA, pa);
    cIasI(pB, pb);
    pa->sign = PLUS;
    pb->sign = PLUS;
    cIasint(xa, 1);
    xa->sign = a->sign;
    cIasint(ya, 0);		/* pa == xa*pA + ya*pB */
    cIasint(xb, 0);
    cIasint(yb, 1);
    yb->sign = b->sign;		/* pb == xb*pA + yb*pB */
    while (Ieven(pa))
    {
	Isr1(pa);
	if (Ieven(xa) && Ieven(ya))
	{
	    Isr1(xa);
	    Isr1(ya);
	}
	else
	{
	    IplasI(xa, pB);
	    Isr1(xa);
	    ImiasI(ya, pA);
	    Isr1(ya);
	}
    }
    while (Ieven(pb))
    {
	Isr1(pb);
	if (Ieven(xb) && Ieven(yb))
	{
	    Isr1(xb);
	    Isr1(yb);
	}
	else
	{
	    IplasI(xb, pB);
	    Isr1(xb);
	    ImiasI(yb, pA);
	    Isr1(yb);
	}
    }
    while (TRUE)
    {
	if (IgtI(pb, pa))
	{
	    swap = pa;
	    pa = pb;
	    pb = swap;
	    swap = xa;
	    xa = xb;
	    xb = swap;
	    swap = ya;
	    ya = yb;
	    yb = swap;
	}
	else
	{
	    ImiasI(pa, pb);
	    ImiasI(xa, xb);
	    ImiasI(ya, yb);
	    if (!Ieq0(pa))
		while (Ieven(pa))
		{
		    Isr1(pa);
		    if (Ieven(xa) && Ieven(ya))
		    {
			Isr1(xa);
			Isr1(ya);
		    }
		    else
		    {
			IplasI(xa, pB);
			Isr1(xa);
			ImiasI(ya, pA);
			Isr1(ya);
		    }
		}
	    else
		break;
	}
    }
    IasIslD(d, pb, toshift);
    IasI(u, xb);
    IasI(v, yb);
    dI(pa);
    dI(pb);
    dI(pA);
    dI(pB);
    dI(xa);
    dI(xb);
    dI(ya);
    dI(yb);
}				/* Ibxgcd */

/***********************************************/

void 
Ireduce(a, b)
    pInteger a, b;
/* Kuerze gcd von a und b */
{
    Integer d, r;
    cI(&d);
    cI(&r);
    Ibgcd(&d, a, b);
    Idiv(a, &r, a, &d);
    Idiv(b, &r, b, &d);
    dI(&d);
    dI(&r);
}