ilib.c

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

/* Integer Version 2.0, TP, 11.02.97,   ilib.c  */ 

#include <iint.h>
#include <math.h> 

void
IasIpowD(x, a, n)
  register pInteger x;
  register const Integer *a;
  DigitType n;
{
  Integer y, h;
  DigitType m;

  if (Ieq0(a))
    Ias0(x);
  else if (Ieq1(a))
    Ias1(x);
  else
  {
    cI(&h);
    cI(&y);
    IasI(&y, a);
    Ias1(x);
    m = n;

    while (m > 0)
    {
      if (m & 1)
	ImuasI(x, &y);
      ImuasI(&y, &y);
      m >>= 1;
    }
    dI(&h);
    dI(&y);
  }
}

/* returns floor(sqrt(DigitType x)) using the double sqrt */
static DigitType 
DigitSqrt (x)
     DigitType x;
{
  return (DigitType) sqrt ((double) x);
}

/* returns floor(sqrt(Integer x)) in y */ 
void
IassqrtI (y, x)
     pInteger y;
     const Integer *x;
{

  /* x < 0  ? ==> error */
  if (Ilt0 (x))
    Ierror ("IassqrtI: negative argument");

  /* special cases: x == 0, x == 1, x single precision */
  if (Ieq0 (x))
    Ias0 (y);
  else if (Ieq1 (x))
    Ias1 (y);
  else if (x->length == 1)
    Iasulong (y, DigitSqrt (x->vec[0]));

  /* the hard case */
  else
    {
      /* xx is a copy of y; q, r are temporary variables */
      Integer xx, q, r;
      int i = x->length;
      DigitType k;

      cIasI (&xx, x);
      cImaxlength (&q, (i + 1) / 2);
      cImaxlength (&r, (i + 1) / 2);
  
      /* first approximation of sqrt(x) in y */
      k = DigitSqrt (xx.vec[i - 1]);
      Iasulong (y, k + 1);
      IslasD (y, (i - 1) * (BitsPerDigit / 2));

      /* newton loop */
      for (;;)
        {
          uIdiv (&q, &r, &xx, y);
          if (IgeI (&q, y))
              break;
          else
            {
              IplasI (y, &q);
              Isr1 (y);
            }
        }
      dI (&xx);
      dI (&r);
      dI (&q);
    }
}

void
Iroot(y, x, n)
    register pInteger y;
    register const Integer *x;
    register DigitType n;
{
    Integer a, b;
    register DigitType nm1;

    if (Ieq0(x))
    {
      Ias0(y);
      return;
    }

    switch (n)
    {
      case 0:
	Iasint(y, 1);
	return;
      case 1:
	IasI(y, x);
      case 2:
	IassqrtI(y, x);
	return;
    }
    cI(&a);
    cI(&b);
    nm1 = n - 1;
    Iasint(y, 1);
    IslasD(y, (Ilog(x) + nm1) / n);
    do
    {
	IasIpowD(&a, y, nm1);
	Idiv(&a, &b, x, &a);
	IasImiI(&a, y, &a);
	IdiasD(&a, n);
	ImiasI(y, &a);
    } while (Igt0(&a));
    IasIpowD(&a, y, n);
    if (IgtI(&a, x))
	Idec(y);
    dI(&a);
    dI(&b);
}