xdouble.cpp

数值算法库for Windows

      z.e = a.e;
      z.normalize();
      return;
   }
   else {
      if (e > a.e+1) {
         z.x = -x;
         z.e = e;
         z.normalize();
         return;
      }

      z.x = a.x*NTL_XD_BOUND_INV - x;
      z.e = e;
      z.normalize();
      return;
   }
}

double log(const xdouble& a)
{
   static double LogBound = log(NTL_XD_BOUND);
   if (a.x <= 0) {
      Error("log(xdouble): argument must be positive");
   }

   return log(a.x) + a.e*LogBound;
}

xdouble xexp(double x)
{
   const double LogBound = log(NTL_XD_BOUND);

   double y = x/LogBound;
   double iy = floor(y+0.5);

   if (iy >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("xdouble: overflow");

   if (iy <= -(1L << (NTL_BITS_PER_LONG-4)))
      Error("xdouble: underflow");


   double fy = y - iy;

   xdouble res;
   res.e = long(iy);
   res.x = exp(fy*LogBound);
   res.normalize();
   return res;
}

/**************  input / output routines **************/


void ComputeLn2(RR&);
void ComputeLn10(RR&);

long ComputeMax10Power()
{
   long old_p = RR::precision();
   RR::SetPrecision(NTL_BITS_PER_LONG);

   RR ln2, ln10;
   ComputeLn2(ln2);
   ComputeLn10(ln10);

   long k = to_long( to_RR(1L << (NTL_BITS_PER_LONG-5)) * ln2 / ln10 );

   RR::SetPrecision(old_p);
   return k;
}


xdouble PowerOf10(const ZZ& e)
{
   static long init = 0;
   static xdouble v10k;
   static long k;

   if (!init) {
      long old_p = RR::precision();
      k = ComputeMax10Power();
      RR::SetPrecision(NTL_DOUBLE_PRECISION);
      v10k = to_xdouble(power(to_RR(10), k)); 
      RR::SetPrecision(old_p);
      init = 1;
   }

   ZZ e1;
   long neg;

   if (e < 0) {
      e1 = -e;
      neg = 1;
   }
   else {
      e1 = e;
      neg = 0;
   }

   long r;
   ZZ q;

   r = DivRem(q, e1, k);

   long old_p = RR::precision();
   RR::SetPrecision(NTL_DOUBLE_PRECISION);
   xdouble x1 = to_xdouble(power(to_RR(10), r));
   RR::SetPrecision(old_p);

   xdouble x2 = power(v10k, q);
   xdouble x3 = x1*x2;

   if (neg) x3 = 1/x3;

   return x3;
}




ostream& operator<<(ostream& s, const xdouble& a)
{
   if (a == 0) {
      s << "0";
      return s;
   }

   long old_p = RR::precision();
   long temp_p = long(log(fabs(log(fabs(a))) + 1.0)/log(2.0)) + 10; 

   RR::SetPrecision(temp_p);

   RR ln2, ln10, log_2_10;
   ComputeLn2(ln2);
   ComputeLn10(ln10);
   log_2_10 = ln10/ln2;
   ZZ log_10_a = to_ZZ(
  (to_RR(a.e)*to_RR(2*NTL_XD_HBOUND_LOG) + log(fabs(a.x))/log(2.0))/log_2_10);

   RR::SetPrecision(old_p);

   xdouble b;
   long neg;

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

   ZZ k = xdouble::OutputPrecision() - log_10_a;

   xdouble c, d;

   c = PowerOf10(to_ZZ(xdouble::OutputPrecision()));
   d = PowerOf10(log_10_a);

   b = b / d;
   b = b * c;

   while (b < c) {
      b = b * 10.0;
      k++;
   }

   while (b >= c) {
      b = b / 10.0;
      k--;
   }

   b = b + 0.5;
   k = -k;

   ZZ B;
   conv(B, b);

   long bp_len = xdouble::OutputPrecision()+10;

   char *bp = NTL_NEW_OP char[bp_len];

   if (!bp) Error("xdouble output: out of memory");

   long len, i;

   len = 0;
   do {
      if (len >= bp_len) Error("xdouble output: buffer overflow");
      bp[len] = DivRem(B, B, 10) + '0';
      len++;
   } while (B > 0);

   for (i = 0; i < len/2; i++) {
      char tmp;
      tmp = bp[i];
      bp[i] = bp[len-1-i];
      bp[len-1-i] = tmp;
   }

   i = len-1;
   while (bp[i] == '0') i--;

   k += (len-1-i);
   len = i+1;

   bp[len] = '\0';

   if (k > 3 || k < -len - 3) {
      // use scientific notation

      if (neg) s << "-";
      s << "0." << bp << "e" << (k + len);
   }
   else {
      long kk = to_long(k);

      if (kk >= 0) {
         if (neg) s << "-";
         s << bp;
         for (i = 0; i < kk; i++) 
            s << "0";
      }
      else if (kk <= -len) {
         if (neg) s << "-";
         s << "0.";
         for (i = 0; i < -len-kk; i++)
            s << "0";
         s << bp;
      }
      else {
         if (neg) s << "-";
         for (i = 0; i < len+kk; i++)
            s << bp[i];
   
         s << ".";
   
         for (i = len+kk; i < len; i++)
            s << bp[i];
      }
   }

   delete [] bp;
   return s;
}

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

   if (!s) Error("bad xdouble 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 xdouble 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 xdouble 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 xdouble input");

   xdouble t1, t2, v;

   if (got1 || got2) {
      conv(t1, a);
      conv(t2, b);
      v = t1/t2;
   }
   else
      v = 1;

   if (sign < 0)
      v = -v;

   if (got_e) {
      if (e_sign < 0) negate(e, e);
      t1 = PowerOf10(e);
      v = v * t1;
   }

   x = v;
   return s;
}


xdouble to_xdouble(const char *s)
{
   long c;
   long sign;
   ZZ a, b;
   long i=0;

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

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

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

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

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

      i++;
      c = s[i];

      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);
            i++;
            c = s[i];
         }
      }
   }

   if (got_dot && !got1 && !got2)  Error("bad xdouble 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;

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

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

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

   xdouble t1, t2, v;

   if (got1 || got2) {
      conv(t1, a);
      conv(t2, b);
      v = t1/t2;
   }
   else
      v = 1;

   if (sign < 0)
      v = -v;

   if (got_e) {
      if (e_sign < 0) negate(e, e);
      t1 = PowerOf10(e);
      v = v * t1;
   }

   return v;
}

NTL_END_IMPL