rr.cpp

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

#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 (NTL_OVERFLOW(p, 1, 0))
      Error("RR: precision too high");

   prec = p;
}

long RR::oprec = 10;

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

   if (NTL_OVERFLOW(p, 1, 0))
      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 >= NTL_OVFBND)
      Error("RR: overflow");

   if (z.e <= -NTL_OVFBND)
      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 >= NTL_OVFBND)
      Error("MakeRR: e too big");

   if (e <= -NTL_OVFBND)
      Error("MakeRR: e too small");

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

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

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

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


static 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 ConvPrec(RR& x, const RR& a, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("ConvPrec: bad precsion");

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

void RoundToPrecision(RR& x, const RR& a, long p)
{
   ConvPrec(x, a, 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 AddPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("AddPrec: bad precsion");

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

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 SubPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("SubPrec: bad precsion");

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

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

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

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

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

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

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


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 MulPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("MulPrec: bad precsion");

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


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 SqrPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("SqrPrec: bad precsion");

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



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 DivPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("DivPrec: bad precsion");

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


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);
   // since k >= 2*prec - bits(a) + 1, T1 has at least 2*prec+1 bits,           
   // thus T1 >= 2^(2*prec)                                                     

   SqrRoot(t.x, T1); // t.x >= 2^prec thus t.x contains the round bit           
   t.e = (a.e - k)/2;
   sqr(T2, t.x);  

   // T1-T2 is the (lower part of the) sticky bit                               
   normalize(z, t, T2 < T1);
}




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

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




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;

   SubPrec(t, a, b, 1);
   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 TruncPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("TruncPrec: bad precsion");

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

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 FloorPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("FloorPrec: bad precsion");

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

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 CeilPrec(RR& x, const RR& a, const RR& b, long p)
{
   if (p < 1 || NTL_OVERFLOW(p, 1, 0))
      Error("CeilPrec: bad precsion");

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

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;
   ConvPrec(t, a, len+a.e);
   xcopy(z, t);
}

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

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


   

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

void ConvPrec(RR& x, const ZZ& 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(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 ConvPrec(RR& x, long 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(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 ConvPrec(RR& x, unsigned long 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(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);