rr.cpp

NTL is a high-performance, portable C++ library providing data structures and algorithms for manipul

   xcopy(z, t);
}

void ConvPrec(RR& x, double a, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("ConvPrec: bad precsion");

   long old_p = RR::prec;
   RR::prec = p;
   conv(x, a);
   RR::prec = old_p;
}


void conv(ZZ& z, const RR& a)
{
   if (a.e >= 0) 
      LeftShift(z, a.x, a.e);
   else {
      long sgn = sign(a.x);
      RightShift(z, a.x, -a.e);
      if (sgn < 0)
         sub(z, z, 1);
   }
}

void CeilToZZ(ZZ& z, const RR& a)
{
   if (a.e >= 0)
      LeftShift(z, a.x, a.e);
   else {
      long sgn = sign(a.x);
      RightShift(z, a.x, -a.e);
      if (sgn > 0)
         add(z, z, 1);
   }
}

void TruncToZZ(ZZ& z, const RR& a)
{
   if (a.e >= 0)
      LeftShift(z, a.x, a.e);
   else 
      RightShift(z, a.x, -a.e);
}


void RoundToZZ(ZZ& z, const RR& a)
{
   if (a.e >= 0) {
      LeftShift(z, a.x, a.e);
      return;
   }

   long len = NumBits(a.x);

   if (-a.e > len) {
      z = 0;
      return;
   }

   if (-a.e == len) {
      if (len == 1)
         z = 0;
      else
         z = sign(a.x);

      return;
   }

   static RR t;

   ConvPrec(t, a, len+a.e);

   LeftShift(z, t.x, t.e);
}


void conv(long& z, const RR& a)
{
   ZZ t;
   conv(t, a);
   conv(z, t);
}

void conv(double& z, const RR& aa)
{
   double x;
   static RR a;

   ConvPrec(a, aa, NTL_DOUBLE_PRECISION);
   // round to NTL_DOUBLE_PRECISION bits to avoid double overflow

   conv(x, a.x);
   z = _ntl_ldexp(x, a.e);
}


void add(RR& z, const RR& a, double b)
{
   static RR B;
   B = b;
   add(z, a, B);
}



void sub(RR& z, const RR& a, double b)
{
   static RR B;
   B = b;
   sub(z, a, B);
}

void sub(RR& z, double a, const RR& b)
{
   static RR A;
   A = a;
   sub(z, A, b);
}



void mul(RR& z, const RR& a, double b)
{
   static RR B;
   B = b;
   mul(z, a, B);
}


void div(RR& z, const RR& a, double b)
{
   static RR B;
   B = b;
   div(z, a, B);
}

void div(RR& z, double a, const RR& b)
{
   static RR A;
   A = a;
   div(z, A, b);
}


void inv(RR& z, const RR& a)
{
   static RR one = to_RR(1);
   div(z, one, a);
}

void InvPrec(RR& x, const RR& a, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("InvPrec: bad precsion");

   long old_p = RR::prec;
   RR::prec = p;
   inv(x, a);
   RR::prec = old_p;
}


long compare(const RR& a, double b)
{
   if (b == 0) return sign(a);

   static RR B;
   B = b;
   return compare(a, B);
}


long operator==(const RR& a, double b) 
{
   if (b == 0) return IsZero(a);
   if (b == 1) return IsOne(a);

   static RR B;
   B = b;
   return a == B;
}


void power(RR& z, const RR& a, long e)
{
   RR b, res;

   long n = NumBits(e);

   long p = RR::precision();
   RR::SetPrecision(p + n + 10);

   xcopy(b, a);

   set(res);
   long i;

   for (i = n-1; i >= 0; i--) {
      sqr(res, res);
      if (bit(e, i))
         mul(res, res, b);
   }

   RR::SetPrecision(p);

   if (e < 0) 
      inv(z, res);
   else
      xcopy(z, res);
}


istream& operator>>(istream& s, RR& x)
{
   long c;
   long cval;
   long sign;
   ZZ a, b;

   if (!s) Error("bad RR input");


   c = s.peek();
   while (IsWhiteSpace(c)) {
      s.get();
      c = s.peek();
   }

   if (c == '-') {
      sign = -1;
      s.get();
      c = s.peek();
   }
   else
      sign = 1;

   long got1 = 0;
   long got_dot = 0;
   long got2 = 0;

   a = 0;
   b = 1;

   cval = CharToIntVal(c);

   if (cval >= 0 && cval <= 9) {
      got1 = 1;

      while (cval >= 0 && cval <= 9) {
         mul(a, a, 10);
         add(a, a, cval);
         s.get();
         c = s.peek();
         cval = CharToIntVal(c);
      }
   }

   if (c == '.') {
      got_dot = 1;

      s.get();
      c = s.peek();
      cval = CharToIntVal(c);

      if (cval >= 0 && cval <= 9) {
         got2 = 1;
   
         while (cval >= 0 && cval <= 9) {
            mul(a, a, 10);
            add(a, a, cval);
            mul(b, b, 10);
            s.get();
            c = s.peek();
            cval = CharToIntVal(c);
         }
      }
   }

   if (got_dot && !got1 && !got2)  Error("bad RR input");

   ZZ e;

   long got_e = 0;
   long e_sign;

   if (c == 'e' || c == 'E') {
      got_e = 1;

      s.get();
      c = s.peek();

      if (c == '-') {
         e_sign = -1;
         s.get();
         c = s.peek();
      }
      else if (c == '+') {
         e_sign = 1;
         s.get();
         c = s.peek();
      }
      else
         e_sign = 1;

      cval = CharToIntVal(c);

      if (cval < 0 || cval > 9) Error("bad RR input");

      e = 0;
      while (cval >= 0 && cval <= 9) {
         mul(e, e, 10);
         add(e, e, cval);
         s.get();
         c = s.peek();
         cval = CharToIntVal(c);
      }
   }

   if (!got1 && !got2 && !got_e) Error("bad RR input");

   RR t1, t2, v;

   long old_p = RR::precision();

   if (got1 || got2) {
      ConvPrec(t1, a, max(NumBits(a), 1));
      ConvPrec(t2, b, NumBits(b));
      if (got_e)
         RR::SetPrecision(old_p + 10);

      div(v, t1, t2);
   }
   else
      set(v);

   if (sign < 0)
      negate(v, v);

   if (got_e) {
      if (e >= NTL_OVFBND) Error("RR input overflow");
      long E;
      conv(E, e);
      if (e_sign < 0) E = -E;
      RR::SetPrecision(old_p + 10);
      power(t1, to_RR(10), E);
      mul(v, v, t1);
      RR::prec = old_p;
   }

   xcopy(x, v);
   return s;
}

void InputPrec(RR& x, istream& s, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("ConvPrec: bad precsion");

   long old_p = RR::prec;
   RR::prec = p;
   s >> x;
   RR::prec = old_p;
}


void conv(RR& z, const xdouble& a)
{
   conv(z, a.mantissa());

   if (a.exponent() >  ((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG))) 
      Error("RR: overlow");

   if (a.exponent() < -((2*NTL_OVFBND)/(2*NTL_XD_HBOUND_LOG))) 
      Error("RR: underflow");

   z.e += a.exponent()*(2*NTL_XD_HBOUND_LOG);

   if (z.e >= NTL_OVFBND)
      Error("RR: overflow");

   if (z.e <= -NTL_OVFBND)
      Error("RR: underflow");
}

void ConvPrec(RR& x, const xdouble& a, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("ConvPrec: bad precsion");

   long old_p = RR::prec;
   RR::prec = p;
   conv(x, a);
   RR::prec = old_p;
}



void conv(xdouble& z, const RR& a)
{
   xdouble x;
   xdouble y;

   conv(x, a.x);
   power2(y, a.e);
   z = x*y;
}
      
void power2(RR& z, long e)
{
   if (e >= NTL_OVFBND)
      Error("RR: overflow");

   if (e <= -NTL_OVFBND)
      Error("RR: underflow");

   set(z.x); 
   z.e = e;
}

void conv(RR& z, const quad_float& a)
{
   static RR hi, lo, res;

   ConvPrec(hi, a.hi, NTL_DOUBLE_PRECISION);
   ConvPrec(lo, a.lo, NTL_DOUBLE_PRECISION);

   add(res, hi, lo);

   z = res;
}

void ConvPrec(RR& x, const quad_float& a, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("ConvPrec: bad precsion");

   long old_p = RR::prec;
   RR::prec = p;
   conv(x, a);
   RR::prec = old_p;
}


void conv(quad_float& z, const RR& a)
{
   long old_p;
   static RR a_hi, a_lo;

   old_p = RR::prec;

   ConvPrec(a_hi, a, NTL_DOUBLE_PRECISION);  // high order bits
   SubPrec(a_lo, a, a_hi, NTL_DOUBLE_PRECISION);  // low order bits

   z = to_quad_float(a_hi.x)*power2_quad_float(a_hi.e) +
       to_quad_float(a_lo.x)*power2_quad_float(a_lo.e);
}

void conv(RR& x, const char *s)
{
   long c;
   long cval;
   long sign;
   ZZ a, b;
   long i = 0;

   if (!s) Error("bad RR input");


   c = s[i];
   while (IsWhiteSpace(c)) {
      i++;
      c = s[i];
   }

   if (c == '-') {
      sign = -1;
      i++;
      c = s[i];
   }
   else
      sign = 1;

   long got1 = 0;
   long got_dot = 0;
   long got2 = 0;

   a = 0;
   b = 1;

   cval = CharToIntVal(c);

   if (cval >= 0 && cval <= 9) {
      got1 = 1;

      while (cval >= 0 && cval <= 9) {
         mul(a, a, 10);
         add(a, a, cval);
         i++;
         c = s[i];
         cval = CharToIntVal(c);
      }
   }

   if (c == '.') {
      got_dot = 1;

      i++;
      c = s[i];
      cval = CharToIntVal(c);

      if (cval >= 0 && cval <= 9) {
         got2 = 1;
   
         while (cval >= 0 && cval <= 9) {
            mul(a, a, 10);
            add(a, a, cval);
            mul(b, b, 10);
            i++;
            c = s[i];
            cval = CharToIntVal(c);
         }
      }
   }

   if (got_dot && !got1 && !got2)  Error("bad RR input");

   ZZ e;

   long got_e = 0;
   long e_sign;

   if (c == 'e' || c == 'E') {
      got_e = 1;

      i++;
      c = s[i];

      if (c == '-') {
         e_sign = -1;
         i++;
         c = s[i];
      }
      else if (c == '+') {
         e_sign = 1;
         i++;
         c = s[i];
      }
      else
         e_sign = 1;


      cval = CharToIntVal(c);

      if (cval < 0 || cval > 9) Error("bad RR input");

      e = 0;
      while (cval >= 0 && cval <= 9) {
         mul(e, e, 10);
         add(e, e, cval);
         i++;
         c = s[i];
         cval = CharToIntVal(c);
      }
   }

   if (!got1 && !got2 && !got_e) Error("bad RR input");

   RR t1, t2, v;

   long old_p = RR::precision();

   if (got1 || got2) {
      ConvPrec(t1, a, max(NumBits(a), 1));
      ConvPrec(t2, b, NumBits(b));
      if (got_e)
         RR::SetPrecision(old_p + 10);

      div(v, t1, t2);
   }
   else
      set(v);

   if (sign < 0)
      negate(v, v);

   if (got_e) {
      if (e >= NTL_OVFBND) Error("RR input overflow");
      long E;
      conv(E, e);
      if (e_sign < 0) E = -E;
      RR::SetPrecision(old_p + 10);
      power(t1, to_RR(10), E);
      mul(v, v, t1);
      RR::prec = old_p;
   }

   xcopy(x, v);
}

void ConvPrec(RR& x, const char *s, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("ConvPrec: bad precsion");

   long old_p = RR::prec;
   RR::prec = p;
   conv(x, s);
   RR::prec = old_p;
}


void ReallyComputeE(RR& res)
{
   long p = RR::precision();
   RR::SetPrecision(p + NumBits(p) + 10);

   RR s, s1, t;

   s = 1;
   t = 1;

   long i;

   for (i = 2; ; i++) {
      add(s1, s, t);
      if (s == s1) break;
      xcopy(s, s1);
      div(t, t, i);
   }

   RR::SetPrecision(p);
   xcopy(res, s);
}

void ComputeE(RR& res)
{
   static long prec = 0;
   static RR e;

   long p = RR::precision();

   if (prec <= p + 10) {
      prec = p + 20;
      RR::SetPrecision(prec);
      ReallyComputeE(e);
      RR::SetPrecision(p);
   }

   xcopy(res, e);
}


void exp(RR& res, const RR& x)
{
   if (x >= NTL_OVFBND || x <= -NTL_OVFBND)
      Error("RR: overflow");

   long p = RR::precision();

   // step 0: write x = n + f, n an integer and |f| <= 1/2
   //    careful -- we want to compute f to > p bits of precision


   RR f, nn;
   RR::SetPrecision(NTL_BITS_PER_LONG);
   round(nn, x);
   RR::SetPrecision(p + 10);
   sub(f, x, nn);
   long n = to_long(nn);

   // step 1: calculate t1 = e^n by repeated squaring

   RR::SetPrecision(p + NumBits(n) + 10);

   RR e;
   ComputeE(e);

   RR::SetPrecision(p + 10);

   RR t1;
   power(t1, e, n); 

   // step 2: calculate t2 = e^f using Taylor series expansion