// --------------------------------------------------------------------------
// This is our standard testfile for checking the LEDA data type bigfloat
// --------------------------------------------------------------------------

#include <LEDA/bigfloat.h>
#include <math.h>
#include <iostream.h>
#include <fstream.h>
#include <LEDA/list.h>
#include <assert.h>
#include <LEDA/floatingpoint.h>


// --------------------------------------------------------------------------
// The following are special test routines only usable under SUNOS
// compilers: CC, g++ and SOLARIS/CC
// --------------------------------------------------------------------------

#ifdef SUNOS

// extern int ieee_flags(const char*, const char*, const char*, char**);
// scheint nur mit g++ unter SUNOS zu klappen

#include <sunmath.h>

void set_ieee_mode(rounding_modes rmode)
{
  char* out;
  switch(rmode)
  {
    case TO_ZERO: ieee_flags("set","direction","tozero",&out); break;
    case TO_P_INF: ieee_flags("set","direction","positive",&out); break;
    case TO_N_INF: ieee_flags("set","direction","negative",&out); break;
    case TO_NEAREST:   ieee_flags("set","direction","nearest",&out); break;
    case EXACT: case TO_INF: break;
  }
}

bigfloat add_debug(const bigfloat& a,const bigfloat& b, rounding_modes rmode) 
{ 
   double a_d = to_double(a);
   double b_d = to_double(b);
   set_ieee_mode(rmode);
   double result_d = a_d + b_d;
   set_ieee_mode(TO_NEAREST);
   bigfloat result_bf = add(a,b,53,rmode);
   if ( (a != bigfloat(a_d)) || (b != bigfloat(b_d)) )
     error_handler(1,"error add_debug: operands not correct"); 
   if (result_bf != bigfloat(result_d))
     error_handler(1,"error add_debug: result incorrect");
   return result_bf;
}
  
bigfloat sub_debug(const bigfloat& a,const bigfloat& b, rounding_modes rmode) 
{ 
   double a_d = to_double(a);
   double b_d = to_double(b);
   set_ieee_mode(rmode);
   double result_d = a_d - b_d;
   set_ieee_mode(TO_NEAREST);
   bigfloat result_bf = sub(a,b,53,rmode);
   if ( (a != bigfloat(a_d)) || (b != bigfloat(b_d)) )
     error_handler(1,"error sub_debug: operands not correct"); 
   if (result_bf != bigfloat(result_d))
     error_handler(1,"error sub_debug: result incorrect");
   return result_bf;
}

bigfloat mul_debug(const bigfloat& a,const bigfloat& b, rounding_modes rmode) 
{
   double a_d = to_double(a);
   double b_d = to_double(b);
   set_ieee_mode(rmode);
   double result_d = a_d * b_d;
   set_ieee_mode(TO_NEAREST);
   bigfloat result_bf = mul(a,b,53,rmode);
   if ( (a != bigfloat(a_d)) || (b != bigfloat(b_d)) )
     error_handler(1,"error mul_debug: operands not correct"); 
   if (result_bf != bigfloat(result_d))
     error_handler(1,"error mul_debug: result incorrect");
   return result_bf;
}

bigfloat div_debug(const bigfloat& a,const bigfloat& b, rounding_modes rmode) 
{ 
   double a_d = to_double(a);
   double b_d = to_double(b);
   set_ieee_mode(rmode);
   double result_d = a_d / b_d;
   set_ieee_mode(TO_NEAREST);
   bigfloat result_bf = div(a,b,53,rmode);
   if ( (a != bigfloat(a_d)) || (b != bigfloat(b_d)) )
     error_handler(1,"error div_debug: operands not correct"); 
   if (result_bf != bigfloat(result_d))
     error_handler(1,"error div_debug: result incorrect");
   return result_bf;
}

#endif




// --------------------------------------------------------------------------
// Now the main check routines follow
// --------------------------------------------------------------------------


main()
{

  // float s7 = sqrt(7);
  // double sqrt7 = s7;
  // cout << "sqrt7: " << to_double(sqrt(bigfloat(7))) << "\n" << flush;


// --------------------------------------------------------------------------
//   first we test special cases
// --------------------------------------------------------------------------

{
  cout << "\nteste Spezialfaelle\n\n";
 
  // Funktionen:
  // isNaN, isSpecial, isZero, ispZero, isnZero, isInf, ispInf, isnInf
  // Arithmetik fuer special values

  assert(ispZero(bigfloat::pZero));
  assert(isnZero(bigfloat::nZero));
  assert(isZero(bigfloat::pZero));
  assert(isZero(bigfloat::nZero));
  assert(isZero(bigfloat::zero));

  assert(ispInf(bigfloat::pInf));
  assert(isnInf(bigfloat::nInf));
  assert(isInf(bigfloat::pInf));
  assert(isInf(bigfloat::nInf));

  assert(isNaN(bigfloat::NaN));


  list<bigfloat> S, T;
  bigfloat x, y;
  double dx, dy;

  S.append(bigfloat::pZero);
  S.append(bigfloat::nZero);
  S.append(bigfloat::pInf);
  S.append(bigfloat::nInf);

  forall(x,S)
  {
     assert (isSpecial(x));

     if (x == bigfloat::zero)
        assert(isZero(x));
     else
        assert(!isZero(x));

     if (x == bigfloat::pInf)
     {
        assert(ispInf(x));
        assert(isInf(x));
     }
     else
        assert(!ispInf(x));

     if (x == bigfloat::nInf)
     {
        assert(isnInf(x));
        assert(isInf(x));
     }
     else
        assert(!isnInf(x));
  }

  y = bigfloat::NaN;
  forall(x,S)
  {
     assert(bigfloat(to_double(x)) == x);
     assert(isNaN(x+y)&&isNaN(x-y)&&isNaN(x*y)&&isNaN(x/y));
     assert(isNaN(y+x)&&isNaN(y-x)&&isNaN(y*x)&&isNaN(y/x));
  }

  T = S;
  forall(x,S) 
     forall(y,T)
  {
     dx = to_double(x);
     dy = to_double(y);

     // Note: Some compilers like sunpro give NaN == NaN for doubles 
     // This is not IEEE-conform, but the fact remains

     assert(   binary_equal(to_double(x+y),dx+dy) 
            || (isNaN(x+y) && is_nan(dx+dy))
           );
     assert(   binary_equal(to_double(x-y),dx-dy) 
            || (isNaN(x-y) && is_nan(dx-dy))
           );
     assert(   binary_equal(to_double(x*y),dx*dy) 
            || (isNaN(x*y) && is_nan(dx*dy))
           );
     assert(   binary_equal(to_double(x/y),dx/dy) 
            || (isNaN(x/y) && is_nan(dx/dy))
           );
  }
 
} // special cases



// --------------------------------------------------------------------------
// Rundungsfunktionen und Konvertierung
// --------------------------------------------------------------------------


{  

  // alter Konvertierungstest

  cout << "teste Konvertieren double-bigfloat\n";

  int i,j;
  double powi = 1; 
  double x, y;
  double small = ldexp(1,-1000);
  double large = ldexp(1, 1000);

  for (i=1;i<=53;i++)
  {
    powi = powi*2;
    x = (powi-1)*small;
    y = -(powi-1)*large;
    for (j=1;j<=52+i;j++)
    { 
      x = x/double(2);
      assert(x == to_double(bigfloat(x)));
      y = 2*y;
      // assert(y == to_double(bigfloat(y)));
    }
  }

}


{
   // neuer Rundungstest
   
   // Funktionen:
   // round, to_double, double-Konstruktor, get_significant_length,
   // pow2, ilog2, sign, set_precision, to_integer
   

   int i, j;
   bigfloat x, y, z;
   double dx, dy;
   float fx;
   long len;
   bool exact;

   cout << "\nteste Runden, Konvertieren, Basisfunktionen\n\n";

   for (i=1; i<= 500; i++)
   {
      dx = sin(i);
      fx = dx;
      x = bigfloat(dx);
      y = round(x,24);
     
      assert(bigfloat(fx) == y);
      assert(dx == to_double(x));

      z = round(y,24,TO_NEAREST,exact,1);
      assert(z == y && !exact);
      z = round(y,24,TO_NEAREST,exact,-1);
      assert(z == y && !exact);
      z = round(y,24,TO_NEAREST,exact,0);
      assert(z == y && exact);

      z = round(y,24,TO_P_INF,exact,1);
      if (z <= y)
      {
        cout << "z = " << to_double(z) << "\n";
        cout << "y = " << to_double(y)<< "\n";
      }
      assert(z > y && !exact);
      z = round(y,24,TO_P_INF,exact,-1);
      assert(z == y && !exact);

      z = round(y,24,TO_N_INF,exact,1);
      assert(z == y && !exact);
      z = round(y,24,TO_N_INF,exact,-1);
      assert(z < y && !exact);

      z = round(y,24,TO_ZERO,exact,1);
      if (sign(z) > 0)
        assert(z == y && !exact);
      else
        assert(z > y && !exact);
      z = round(y,24,TO_ZERO,exact,-1);
      if (sign(z) > 0)
        assert(z < y && !exact);
      else
        assert(z == y && !exact);

      z = round(y,24,TO_INF,exact,1);
      if (sign(z) > 0)
        assert(z > y && !exact);
      else
        assert(z == y && !exact);
      z = round(y,24,TO_INF,exact,-1);
      if (sign(z) > 0)
        assert(z == y && !exact);
      else
        assert(z < y && !exact);
      
      assert(abs(x) == sign(x)*x);
      len = x.get_significant_length()-x.get_significant().zeros();
      bigfloat::set_precision(len-1);
      assert(abs(x) != sign(x)*x);
      bigfloat::set_precision(53);
      
      assert(sign(x) == sign(x));
      assert(i == ilog2(pow2(integer(i))).tolong());

      assert(floor(1198357*dx) == floor(1198357*x));
      assert(ceil(1198357*dx) == ceil(1198357*x));
      assert(floor(1198357*dx) <= to_integer(1198357*x));
      assert(to_integer(1198357*x) <= ceil(1198357*x));
   }



   long logx;
   cout << "i = ";
   for (i=1; i<= 50; i++)
   {

   cout << i << " " << flush;   
  
   for (j=1; j<= 50; j++)
   {
      dx = sin(i);
      x = bigfloat(dx);
      y = round(x,j);

      assert(y.get_significant_length() <= j);

      assert(y == bigfloat(to_double(y)));
      assert(y == bigfloat(y.get_significant(),y.get_exponent()));

      assert(ilog2(x) == x.get_significant_length()+x.get_exponent());
      logx = ilog2(x).tolong();

      assert(y == bigfloat(to_double(y))); 

      y = round(x,j,TO_ZERO);
      assert(y.get_significant_length() <= j);
      dy = to_double(y);
      assert (fabs(dy) <= fabs(dx));
      assert (fabs(dx-dy) <= pow(2,logx-j));

      y = round(x,j,TO_INF);
      assert(y.get_significant_length() <= j);
      dy = to_double(y);
      assert (fabs(dy) >= fabs(dx));
      assert (fabs(dx-dy) <= pow(2,logx-j));

      y = round(x,j,TO_N_INF);
      assert(y.get_significant_length() <= j);
      dy = to_double(y);
      assert (dy <= dx);
      assert (fabs(dx-dy) <= pow(2,logx-j));

      y = round(x,j,TO_P_INF);
      assert(y.get_significant_length() <= j);
      dy = to_double(y);
      assert (dy >= dx);
      assert (fabs(dx-dy) <= pow(2,logx-j));

      y = round(x,j,TO_NEAREST);
      assert(y.get_significant_length() <= j);
      dy = to_double(y);
      assert (fabs(dx-dy) <= pow(2,logx-j-1));
       
   } // for j
 
   } // for i

   bigfloat::set_precision(53);

} // Runden und Konvertieren


{
  // Holger Kappels Test mit to_integer und floor

  int i;
  bigfloat x;
  double x_dbl;

  for(i=-100; i<100; i++)
  {
    x = i;
    x = x/bigfloat(50);
    x_dbl = to_double(x);

    assert(2*abs(to_integer(x)-x) <= 1);
    assert (floor(x) <= x && x-floor(x) < 1);
    assert(ceil(x) >= x  && ceil(x)-x  < 1);
    assert(x <= bigfloat(i+1)/bigfloat(50));
    assert(x > bigfloat(i-1)/bigfloat(50));
  }
}



// -------------------------------------------------------------------------
//  Arithmetik
// -------------------------------------------------------------------------

cout << "\n\nteste Arithmetik\n";

{

 // alter Arithmetik-Test

 double x, y;
 int exp;
 bigfloat x_bf, y_bf;
 int i, j;

 cout << "\ni = ";
 for (i=1;i<=25;i++)
 { 
   cout << i << " " << flush;
   for (j=1;j<=i;j++)
     for (exp=0;exp<=53 + j - i;exp++)
   {
     x = (pow(2,i)-1)*pow(2,exp);
     y = pow(2,j)-1;

     // x = 11  ...  1 * 2^exp
     //     [i times ]

     // y = 11  ...  1
     //     [j times ]

     x_bf = x; y_bf = y;

     assert(x+y == to_double(add(x_bf,y_bf,53,TO_NEAREST)));
     assert(y-x == to_double(sub(y_bf,x_bf,53,TO_NEAREST)));
     if (exp == 0) // here exp not important
        assert(x*y == to_double(mul(x_bf,y_bf,53,TO_NEAREST)));
     if (exp == 0) // here exp not important
        assert(x/y == to_double(div(x_bf,y_bf,53,TO_NEAREST)));
   } // for j, exp

 } // for i

}

cout << "\ni = ";


// Funktionen: add, sub, mul, div in allen rounding modes
// Vergleichsoperatoren, set_rounding_mode

{  // neuer Arithmetik-Test

   bool exact;

   int i, j;
   bigfloat x, y;
   double dx, dy;
   long logx;

   bigfloat result;
   bigfloat exact_result;
   bigfloat rounding;
   for (i=1; i<= 50; i++)
   {
   cout << i << " " << flush;
      
   for (j=1; j<= 50; j++)
   {
      dx = sin(i);
      x = bigfloat(dx);

      dy = cos(j);
      y = bigfloat(dy);

 // conversion:
      assert(to_double(x) == dx);
      assert(to_double(y) == dy);

 // comparison:

      assert((x > y)  == (dx > dy));
      assert((x < y)  == (dx < dy));
      assert((x >= y) == (dx >= dy));
      assert((x <= y) == (dx <= dy));

 // addition:

      exact_result = add(x,y,53,EXACT);
      assert(to_double(exact_result) == dx+dy);
      assert(round(exact_result,53,TO_NEAREST) == dx+dy);
      result = add(x,y,53,TO_NEAREST,exact);
      assert(result == dx+dy);
      assert(to_double(result) == dx+dy);
      if (result == exact_result)
         assert(exact);
      else
         assert(!exact);
      logx = ilog2(result).tolong();

      rounding = add(x,y,53,TO_P_INF,exact);
      assert(rounding >= result);
      assert(abs(rounding-exact_result) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);

      rounding = add(x,y,53,TO_N_INF);
      assert(rounding <= result);
      assert(abs(exact_result-rounding) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);

      rounding = add(x,y,53,TO_ZERO);
      assert(abs(rounding) <= abs(result));
      assert(abs(rounding-exact_result) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);

      rounding = add(x,y,53,TO_INF);
      assert(abs(rounding) >= abs(result));
      assert(abs(rounding-exact_result) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);

 // subtraction:

      result = sub(x,y,53,TO_NEAREST);
      assert(to_double(result) == dx-dy);
      logx = ilog2(result).tolong();

      rounding = sub(x,y,53,TO_P_INF);
      assert(rounding >= result);
      assert(abs(rounding-result) <= pow(2,logx-53));

      rounding = sub(x,y,53,TO_N_INF);
      assert(rounding <= result);
      assert(abs(result-rounding) <= pow(2,logx-53));

      rounding = sub(x,y,53,TO_ZERO);
      assert(abs(rounding) <= abs(result));
      assert(abs(rounding-result) <= pow(2,logx-53));

      rounding = sub(x,y,53,TO_INF);
      assert(abs(rounding) >= abs(result));
      assert(abs(rounding-result) <= pow(2,logx-53));

 // multiplication:

      exact_result = mul(x,y,53,EXACT);
      result = mul(x,y,53,TO_NEAREST,exact);
      assert(to_double(result) == dx*dy);
      if (result == exact_result)
         assert(exact);
      else
         assert(!exact);
      logx = ilog2(result).tolong();

      rounding = mul(x,y,53,TO_P_INF,exact);
      assert(rounding >= result);
      assert(abs(rounding-exact_result) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);

      rounding = mul(x,y,53,TO_N_INF,exact);
      assert(rounding <= result);
      assert(abs(exact_result-rounding) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);

      rounding = mul(x,y,53,TO_ZERO,exact);
      assert(abs(rounding) <= abs(result));
      assert(abs(rounding-exact_result) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);

      rounding = mul(x,y,53,TO_INF,exact);
      assert(abs(rounding) >= abs(result));
      assert(abs(rounding-exact_result) <= pow(2,logx-53));
      if (rounding == exact_result)
         assert(exact);
      else
         assert(!exact);



 // division:

      result = div(x,y,53,TO_NEAREST);
      assert(to_double(result) == dx/dy);
      logx = ilog2(result).tolong();

      rounding = div(x,y,53,TO_P_INF);
      assert(rounding >= result);
      assert(abs(rounding-result) <= pow(2,logx-53));

      rounding = div(x,y,53,TO_N_INF);
      assert(rounding <= result);
      assert(abs(result-rounding) <= pow(2,logx-53));

      rounding = div(x,y,53,TO_ZERO);
      assert(abs(rounding) <= abs(result));
      assert(abs(rounding-result) <= pow(2,logx-53));

      rounding = div(x,y,53,TO_INF);
      assert(abs(rounding) >= abs(result));
      assert(abs(rounding-result) <= pow(2,logx-53));

   } // for j

   } // for i

   result = div(bigfloat(4),bigfloat(2),2,TO_NEAREST,exact);
   assert(result == 2 && exact);
   result = div(bigfloat(4),bigfloat(3),2,TO_NEAREST,exact);
   assert(!exact);

   result = div(bigfloat(4),bigfloat(2),2,TO_P_INF,exact);
   assert(result == 2 && exact);
   result = div(bigfloat(4),bigfloat(3),2,TO_P_INF,exact);
   assert(!exact);

   result = div(bigfloat(4),bigfloat(2),2,TO_N_INF,exact);
   assert(result == 2 && exact);
   result = div(bigfloat(4),bigfloat(3),2,TO_N_INF,exact);
   assert(!exact);
  
   result = div(bigfloat(4),bigfloat(2),2,TO_INF,exact);
   assert(result == 2 && exact);
   result = div(bigfloat(4),bigfloat(3),2,TO_INF,exact);
   assert(!exact);

   result = div(bigfloat(4),bigfloat(2),2,TO_ZERO,exact);
   assert(result == 2 && exact);
   result = div(bigfloat(4),bigfloat(3),2,TO_ZERO,exact);
   assert(!exact);

   result = sub(x,x,2,TO_NEAREST,exact);
   assert(result == 0 && EXACT);
   result = sub(x,x,2,TO_N_INF,exact);
   assert(result == 0 && EXACT);
   result = sub(x,x,2,TO_P_INF,exact);
   assert(result == 0 && EXACT);
   result = sub(x,x,2,TO_ZERO,exact);
   assert(result == 0 && EXACT);
   result = sub(x,x,2,TO_INF,exact);
   assert(result == 0 && EXACT);

   x = 1;

   bigfloat::set_rounding_mode(TO_P_INF);
   y = x + pow2(-100);
   assert(y > x);

   bigfloat::set_rounding_mode(TO_INF);
   y = x + pow2(-100);
   assert(y > x);

   bigfloat::set_rounding_mode(TO_N_INF);
   y = x - pow2(-100);
   assert(y < x);

   bigfloat::set_rounding_mode(TO_ZERO);
   y = x - pow2(-100);
   assert(y < x);

   bigfloat::set_rounding_mode(EXACT);
   y = x + pow2(-100);
   assert(y > x);

   bigfloat::set_rounding_mode(TO_NEAREST);
   
   cout << "\n";

} // Arithmetik


// ---------------------------------------------------------------------------
// test sqrt with different precisions and rounding modes
// ---------------------------------------------------------------------------

{

 // alter und neuer Test gemischt

 // Funktion: sqrt

 cout << "\nteste Wurzelfunktion\n\n";

 int prec;
 bigfloat seven = 7, sqrt7;
 bigfloat sqrt_low,sqrt_high;
 bigfloat sum;
 integer exponent, sig;

 cout << "prec = ";

 for (prec=2;prec<=10000;prec*=2)
 {
   cout << prec << " " << flush;

   sqrt7 = sqrt(seven,prec,TO_NEAREST);
   exponent = sqrt7.get_exponent();
   sig = sqrt7.get_significant()-1;
   sqrt_low  = bigfloat(sig,exponent);
   sig = sqrt7.get_significant()+1;
   sqrt_high = bigfloat(sig,exponent);

   // first we test whether sqrt_low is less than sqrt(7) and whether
   // sqrt_high is more than sqrt(7)
 
   assert(mul(sqrt_low ,sqrt_low ,53,EXACT) < seven);
   assert(mul(sqrt_high,sqrt_high,53,EXACT) > seven);

   // now we test whether sqrt7 really is nearer to sqrt(7) than 
   // sqrt_low and sqrt_high

   if  (mul(sqrt7,sqrt7,53,EXACT) >= seven)
   {
     // compare sqrt7 - sqrt(7) and sqrt(7) - sqrt_low
     sum = add(sqrt7,sqrt_low,1,EXACT);
     assert(mul(sum,sum,53,EXACT) <= bigfloat(28));
   }
   else
   {
     // compare sqrt(7) - sqrt7 and sqrt_high - sqrt(7)
     sum = add(sqrt7,sqrt_high,1,EXACT);
     assert(mul(sum,sum,53,EXACT) >= bigfloat(28));
   }

 } // for 

 bigfloat root, next, square;
 int m;
 for (prec = 2; prec <= 100; prec+=10)
 {
     m = prec;

     root = sqrt(m,prec,TO_ZERO);
     square = mul(root,root,1,EXACT);
     assert(square <= m);
     next = bigfloat(root.get_significant()+1,root.get_exponent());
     square = mul(next,next,1,EXACT);
     assert(square > m);

     root = sqrt(m,prec,TO_N_INF);
     square = mul(root,root,1,EXACT);
     assert(square <= m);
     next = bigfloat(root.get_significant()+1,root.get_exponent());
     square = mul(next,next,1,EXACT);
     assert(square > m);

     root = sqrt(m,prec,TO_INF);
     square = mul(root,root,1,EXACT);
     assert(square >= m);
     next = bigfloat(root.get_significant()-1,root.get_exponent());
     square = mul(next,next,1,EXACT);
     assert(square < m);

     root = sqrt(m,prec,TO_P_INF);
     square = mul(root,root,1,EXACT);
     assert(square >= m);
     next = bigfloat(root.get_significant()-1,root.get_exponent());
     square = mul(next,next,1,EXACT);
     assert(square < m);
 }

 bool exact;
 root = sqrt(bigfloat(4),2,TO_NEAREST,exact);
 assert(root == 2);
 assert(exact);
 root = sqrt(bigfloat(5),2,TO_NEAREST,exact);
 assert(!exact);

 root = sqrt(bigfloat(4),2,TO_ZERO,exact);
 assert(root == 2 & exact);
 root = sqrt(bigfloat(5),2,TO_ZERO,exact);
 assert(!exact);

 root = sqrt(bigfloat(4),2,TO_N_INF,exact);
 assert(root == 2 & exact);
 root = sqrt(bigfloat(5),2,TO_N_INF,exact);
 assert(!exact);

 root = sqrt(bigfloat(4),2,TO_P_INF,exact);
 assert(root == 2 & exact);
 root = sqrt(bigfloat(5),2,TO_P_INF,exact);
 assert(!exact);

 root = sqrt(bigfloat(4),2,TO_INF,exact);
 assert(root == 2 & exact);
 root = sqrt(bigfloat(5),2,TO_INF,exact);
 assert(!exact);

} //sqrt


// ---------------------------------------------------------------------------
// input/output
// ---------------------------------------------------------------------------

{
   cout << "\n\nteste Ein- und Ausgabe\n";
   
   bigfloat x, y;

   for (int i=-100;i<=100;i+=5)
   {
      ofstream ostrm("dummy",ios::out);
      x = sin(i);
      ostrm << x << flush;

      ifstream istrm("dummy",ios::in);
      istrm >> y;

      assert(abs(x-y) <= pow2(ilog2(x)-40));
   }

   bigfloat one_half = 1/2.0;
   cout << "1/2: " << one_half << "\n";
  
}

#ifdef SUNOS

// -------------------------------------------------------------------------
// tests with various operations and all rounding modes of IEEE, 
// in precision 53; we use the taylor series to approximate sin(x)
// NOTE: works only with /KM/local/bin/g++ under SUNOS and switch 
// SUNOS must be defined.
// -------------------------------------------------------------------------

{

  cout << "tests with various operations and all rounding modes of ";
  cout << "IEEE,in precision 53;\n we use the taylor series to";
  cout << "approximate sin(x)\n";

  bigfloat::set_precision(53);

  bigfloat s, ds, z, x_sq, sign = -1;
  int k ;

  int rmode;

  for (j=1;j<=50;j++) 
  {
    cout << j << " " << flush;
    x_sq = j*j;
 
    for (rmode = TO_NEAREST; rmode <= TO_N_INF; rmode++)
    {
       bigfloat::set_rounding_mode(TO_NEAREST);

       s = j; ds = j; sign = -1; 
       i = 1;
       while (ds >= ldexp(1,-60))
       {
         ds = mul_debug(ds,x_sq,(rounding_modes) rmode);
         k  = 2*i*(2*i+1);
         ds = div(ds,k,(rounding_modes) rmode);
         s  = add(s,sign*ds, (rounding_modes) rmode);
         sign = -sign;
        i++;
       }
    }
  }

  cout << "\n";


}

#endif

// ----------------------------------------------------------------------------
//  fertig
// ----------------------------------------------------------------------------

cout << "\nAlle Tests absolviert\n\n";

return 0;

} // main
