rr.cpp

数值算法库for Windows

#include <NTL/RR.h>


#include <NTL/new.h>

NTL_START_IMPL


long RR::prec = 150;

void RR::SetPrecision(long p)
{
   if (p < 53)
      p = 53;

   if (p >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: precision too high");

   prec = p;
}

long RR::oprec = 10;

void RR::SetOutputPrecision(long p)
{
   if (p < 1)
      p = 1;

   if (p >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: output precision too high");

   oprec = p;
}



static
void normalize1(RR& z, const ZZ& y_x, long y_e, long prec, long residual)
{
   long len = NumBits(y_x);

   if (len > prec) {
      long correction = ZZ_RoundCorrection(y_x, len - prec, residual);

      RightShift(z.x, y_x, len - prec);

      if (correction) 
         add(z.x, z.x, correction);

      z.e = y_e + len - prec;
   }
   else if (len == 0) {
      clear(z.x);
      z.e = 0;
   }
   else {
      z.x = y_x;
      z.e = y_e;
   }

   if (!IsOdd(z.x))
      z.e += MakeOdd(z.x);

   if (z.e >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: overflow");

   if (z.e <= -(1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: underflow");
}

void normalize(RR& z, const RR& y, long residual = 0)
{
   normalize1(z, y.x, y.e, RR::prec, residual);
}

void MakeRR(RR& z, const ZZ& a,  long e)
{
   if (e >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("MakeRR: e too big");

   if (e <= -(1L << (NTL_BITS_PER_LONG-4)))
      Error("MakeRR: e too small");

   normalize1(z, a, e, RR::prec, 0);
}

void random(RR& z)
{
   static RR t;
   RandomBits(t.x, RR::prec); 
   t.e = -RR::prec;
   normalize(z, t);
}


inline void xcopy(RR& x, const RR& a)
   { normalize(x, a); }

// xcopy emulates old assignment semantics...
// many routines here implicitly assume assignment normalizes,
// but that is no longer the case as of v3.0.

void RoundToPrecision(RR& x, const RR& a, long p)
{
   if (p < 1 || p >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("RoundToPrecision: bad precsion");

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

void conv(RR& x, const RR& a)
{
   normalize(x, a);
}


long IsZero(const RR& a)
{
   return IsZero(a.x);
}

long IsOne(const RR& a)
{
   return a.e == 0 && IsOne(a.x);
}

long sign(const RR& a)
{
   return sign(a.x);
}

void clear(RR& z)
{
   z.e = 0;
   clear(z.x);
}

void set(RR& z)
{
   z.e = 0;
   set(z.x);
}


void add(RR& z, const RR& a, const RR& b)
{
   static RR t;

   if (IsZero(a.x)) {
      xcopy(z, b);
      return;
   }

   if (IsZero(b.x)) {
      xcopy(z, a);
      return;
   }

   if (a.e > b.e) {
      if (a.e-b.e - max(RR::prec-NumBits(a.x),0) >= NumBits(b.x) + 2)
         normalize(z, a, sign(b));
      else {
         LeftShift(t.x, a.x, a.e-b.e);
         add(t.x, t.x, b.x);
         t.e = b.e;
         normalize(z, t);
      }
   }
   else if (a.e < b.e) {
      if (b.e-a.e - max(RR::prec-NumBits(b.x),0) >= NumBits(a.x) + 2)
         normalize(z, b, sign(a));
      else {
         LeftShift(t.x, b.x, b.e-a.e);
         add(t.x, t.x, a.x);
         t.e = a.e;
         normalize(z, t);
      }
   }
   else {
      add(t.x, a.x, b.x);
      t.e = a.e;
      normalize(z, t);
   }
}

void sub(RR& z, const RR& a, const RR& b)
{
   static RR t;

   if (IsZero(a.x)) {
      negate(z, b);
      return;
   }

   if (IsZero(b.x)) {
      xcopy(z, a);
      return;
   }

   if (a.e > b.e) {
      if (a.e-b.e - max(RR::prec-NumBits(a.x),0) >= NumBits(b.x) + 2)
         normalize(z, a, -sign(b));
      else {
         LeftShift(t.x, a.x, a.e-b.e);
         sub(t.x, t.x, b.x);
         t.e = b.e;
         xcopy(z, t);
      }
   }
   else if (a.e < b.e) {
      if (b.e-a.e - max(RR::prec-NumBits(b.x),0) >= NumBits(a.x) + 2) {
         normalize(z, b, -sign(a));
         negate(z.x, z.x);
      }
      else {
         LeftShift(t.x, b.x, b.e-a.e);
         sub(t.x, a.x, t.x);
         t.e = a.e;
         xcopy(z, t);
      }
   }
   else {
      sub(t.x, a.x, b.x);
      t.e = a.e;
      normalize(z, t);
   }
}

void negate(RR& z, const RR& a)
{
   xcopy(z, a);
   negate(z.x, z.x);
}

void abs(RR& z, const RR& a)
{
   xcopy(z, a);
   abs(z.x, z.x);
}


void mul(RR& z, const RR& a, const RR& b)
{
   static RR t;

   mul(t.x, a.x, b.x);
   t.e = a.e + b.e;
   xcopy(z, t);
}

void sqr(RR& z, const RR& a)
{
   static RR t;

   sqr(t.x, a.x);
   t.e = a.e + a.e;
   xcopy(z, t);
}

void div(RR& z, const RR& a, const RR& b)
{
   if (IsZero(b))
      Error("RR: division by zero");

   if (IsZero(a)) {
      clear(z);
      return;
   }

   long la = NumBits(a.x);
   long lb = NumBits(b.x);

   long neg = (sign(a) != sign(b));

   long k = RR::prec - la + lb + 1;
   if (k < 0) k = 0;

   static RR t;
   static ZZ A, B, R;

   abs(A, a.x);
   LeftShift(A, A, k);

   abs(B, b.x);
   DivRem(t.x, R, A, B);

   t.e = a.e - b.e - k;

   normalize(z, t, !IsZero(R));

   if (neg)
      negate(z.x, z.x);
}

void SqrRoot(RR& z, const RR& a)
{
   if (sign(a) < 0)
      Error("RR: attempt to take square root of negative number");

   if (IsZero(a)) {
      clear(z);
      return;
   }

   RR t;
   ZZ T1, T2;
   long k;

   k = 2*RR::prec - NumBits(a.x) + 1;

   if (k < 0) k = 0;

   if ((a.e - k) & 1) k++;

   LeftShift(T1, a.x, k);
   SqrRoot(t.x, T1);
   t.e = (a.e - k)/2;
   sqr(T2, T1);

   normalize(z, t, T2 < T1);
}
   




void swap(RR& a, RR& b)
{
   swap(a.x, b.x);
   swap(a.e, b.e);
}

long compare(const RR& a, const RR& b)
{
   static RR t;

   sub(t, a, b);
   return sign(t);
}



long operator==(const RR& a, const RR& b) 
{
   return a.e == b.e && a.x == b.x;
}


void trunc(RR& z, const RR& a)
{
   static RR t;

   if (a.e >= 0) 
      xcopy(z, a);
   else {
      RightShift(t.x, a.x, -a.e);
      t.e = 0;
      xcopy(z, t);
   }
}

void floor(RR& z, const RR& a)
{
   static RR t;

   if (a.e >= 0) 
      xcopy(z, a);
   else {
      RightShift(t.x, a.x, -a.e);
      if (sign(a.x) < 0)
         add(t.x, t.x, -1);
      t.e = 0;
      xcopy(z, t);
   }
}

void ceil(RR& z, const RR& a)
{
   static RR t;

   if (a.e >= 0)
      xcopy(z, a);
   else {
      RightShift(t.x, a.x, -a.e);
      if (sign(a.x) > 0)
         add(t.x, t.x, 1);
      t.e = 0;
      xcopy(z, t);
   }
}

void round(RR& z, const RR& a)
{
   if (a.e >= 0) {
      xcopy(z, a);
      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;
   RoundToPrecision(t, a, len+a.e);
   xcopy(z, t);
}


   

void conv(RR& z, const ZZ& a)
{
   normalize1(z, a, 0, RR::prec, 0);
}

void conv(RR& z, long a)
{
   if (a == 0) {
      clear(z);
      return;
   }

   if (a == 1) {
      set(z);
      return;
   }

   static ZZ t;
   t = a;
   conv(z, t);
}

void conv(RR& z, unsigned long a)
{
   if (a == 0) {
      clear(z);
      return;
   }

   if (a == 1) {
      set(z);
      return;
   }

   static ZZ t;
   conv(t, a);
   conv(z, t);
}


void conv(RR& z, double a)
{
   if (a == 0) {
      clear(z);
      return;
   }

   if (a == 1) {
      set(z);
      return;
   }

   if (!IsFinite(&a))
      Error("RR: conversion of a non-finite double");

   int e;
   double f;
   static RR t;

   f = frexp(a, &e);

   f = f * NTL_FDOUBLE_PRECISION;
   f = f * 4;

   conv(t.x, f);
   t.e = e - (NTL_DOUBLE_PRECISION + 1);

   xcopy(z, t);
}


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;

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

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


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

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

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

   e = a.e;
   if (e != a.e) Error("RR: overflow in conversion to double");
   conv(x, a.x);

   z = ldexp(x, 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);
}


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 neg;

   long n = NumBits(e);

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

   xcopy(b, a);
   if (e < 0) {
      e = -e;
      neg = 1;
   }
   else
      neg = 0;

   set(res);
   long i;

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

   if (neg) 
      inv(z, res);
   else
      xcopy(z, res);

   RR::SetPrecision(p);
}


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

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


   c = s.peek();
   while (c == ' ' || c == '\n' || c == '\t') {
      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;

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

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

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

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

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

   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;

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

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

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

   RR t1, t2, v;

   long old_p = RR::precision();

   if (got1 || got2) {
      RR::SetPrecision(max(NumBits(a), NumBits(b)));
      conv(t1, a);
      conv(t2, b);
      if (got_e)
         RR::SetPrecision(old_p + 10);
      else
         RR::prec = old_p;
      div(v, t1, t2);
   }
   else
      set(v);

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

   if (got_e) {
      if (e >= (1L << (NTL_BITS_PER_LONG-4))) 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 conv(RR& z, const xdouble& a)
{
   conv(z, a.mantissa());

   // the correctness of the following depends critically on
   // the invariant:
   //    abs(z.e) < (1L << (NTL_BITS_PER_LONG-4)) 

   if (a.exponent() >  ((1L << (NTL_BITS_PER_LONG-3))/(2*NTL_XD_HBOUND_LOG))) 
      Error("RR: overlow");

   if (a.exponent() < -((1L << (NTL_BITS_PER_LONG-3))/(2*NTL_XD_HBOUND_LOG))) 
      Error("RR: underflow");

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

   if (z.e >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: overflow");

   if (z.e <= -(1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: underflow");
}


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 >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: overflow");

   if (e <= -(1L << (NTL_BITS_PER_LONG-4)))
      Error("RR: underflow");

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

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

   long old_p = RR::prec;
   RR::prec = NTL_DOUBLE_PRECISION;

   conv(hi, a.hi);
   conv(lo, a.lo);

   RR::prec = old_p;

   add(res, hi, lo);

   z = res;
}


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