#include <LEDA/real.h>
#include <math.h>
#include <LEDA/array.h>
#include <assert.h>

int main()
{
  float time;
  float overall_time;

// ------------------------------------------------------------------------
// teste real::sign an Beispiel aus Voronoidiagrammrechnung
// ------------------------------------------------------------------------

{

cout << "Teste real::sign an Beispiel aus Voronoidiagrammberechnung\n";

real a1= 15,a2= 21,b1= -19, b2= -a2*b1/a1;
real F,G,SQR;

assert((a1*b2+a2*b1).sign() == 0);

SQR= sqrt( (a1*a1+b1*b1)*(a2*a2+b2*b2) );
F= real(2)*(SQR-a1*a2+b1*b2);
G= real(2)*(SQR+a1*a2-b1*b2);

assert ((G.sign() == 1) && (F.sign() == 0));

}


// -------------------------------------------------------------------------
// teste real::sign an Beispiel der Exponentialfunktion (berechne e)
// teste au"serdem, ob Exaktheit erkannt wird.
// -------------------------------------------------------------------------

{

cout << "Teste real::sign, resp real::operator== am Beispiel der e-Funktion mit Taylorreihe\n";

int k = 20;
int i,j;
array<real> euler1(0,k+1), euler2(0,k+1);
real faculty=1;
real fraction;
real y = 1;
real one = 1;

i = 1;
cout << "i = ";
overall_time = used_time();
for (i=1; i <= k; i++)
{
  cout << i << " " << flush;
  faculty *= real(i);
  euler1[i] = euler2[i] = one;
  j = i;
  fraction = 1;
  while (j > 0)
  {
    fraction *= real(j);
    euler1[i] = euler1[i]*real(y) + fraction;
    euler2[i] = euler2[i]*real(y)/real(j) + real(1);
    j--;
  }
  euler1[i] = euler1[i] / faculty;
  time = used_time();
  assert (euler1[i] == euler2[i]);
  if (i == 20)
  {
    cout << "\nZeit fuer Vergleich gleicher Werte: " << used_time() - time;
    cout << "sec \n";
  }
  time = used_time();
  if (i > 1)
    assert (euler1[i] != euler2[i-1]);
  if (i == 20)
  {
    cout << "Zeit fuer Vergleich ungleicher Werte: " << used_time() - time;
    cout << "sec \n";
  }
}
cout << "overall time: " << used_time() - overall_time << "sec \n";


cout << "\nTeste, ob Exaktheit automatisch erkannt wird am Beispiel der e-Funktion mit Taylorreihe\n";

euler1[20] = one;
j = k;
fraction = 1;
time = used_time();
while (j > 0)
{
  fraction *= real(j);
  euler1[20] = euler1[20]*real(y) + fraction;
  j--;
}
cout << "Zeit fuer Taylorauswertung: " << used_time() - time << "sec \n";

time = used_time();
euler1[20].improve_approximation_to(1000);
cout << "Zeit fuer Euler-improve: " << used_time() - time << "sec \n";

assert (isZero(euler1[20].get_bigfloat_error())); 

}


// ---------------------------------------------------------------------------
// teste real::sign an Beispiel der Wurzelberechnung mit Newton
// ---------------------------------------------------------------------------

{

cout << "Teste real::sign, real::improve, real::compute_with_precision, real::compute_with_precision am Beispiel der Wurzelberechnung mit Newton\n";

int prec = 53;

real result;
real it; // iterator in Newton iteration
real& ne = result; // next iteration value in Newton iteration
real diff;

bigfloat sqrt2;
  
cout << "prec = ";
while (prec < 1000)
{
  cout << prec << " " << flush;
  ne = sqrt(2);
  time = used_time();
  do {  it=ne; ne = (it*it + 2)/(2*it); }
  while (abs(ne-it) >= pow2(-prec-1));
  cout << "Zeit fuer Newtoniteration: " << used_time() - time << "sec \n";

  sqrt2 = sqrt(2,prec+1);

  assert(abs(result - sqrt2) < pow2(-prec));

  result.compute_with_precision(53);
  
  time = used_time();
  result.improve_approximation_to(prec+1);
  cout << "Zeit fuer improve: " << used_time() - time << "sec \n";
  assert(abs(to_bigfloat(result)-sqrt2) < pow2(-prec));
  
  result.compute_with_precision(53);
  time = used_time();
  result.compute_with_precision(prec+8);
  cout << "Zeit fuer compute_with_precision: " << used_time() - time << "sec \n";
  assert(abs(to_bigfloat(result)-sqrt2) < pow2(-prec));
 
  result.compute_with_precision(53);
  time = used_time();
  result.compute_with_precision(prec+2);
  cout << "Zeit fuer compute_with_precision: " << used_time() - time << "sec \n";
  assert(abs(to_bigfloat(result)-sqrt2) < pow2(-prec));

  prec*=2;
}
cout << "\n";

}


// --------------------------------------------------------------------------
// interne Funktion detect_degeneracies getestet
// --------------------------------------------------------------------------

{

cout << "teste interne Funktion detect_degeneracies\n";
int m;
real sum=0;
time=used_time();
for (m=1;m<=50;m++) 
{ sum = sum + (2*sqrt(real(2))-sqrt(real(8))); 
} 
cout << "Zeit fuer Summenberechnung: " << used_time() - time << "sec \n";
time=used_time();
assert(sum.sign() == 0);
cout << "Zeit fuer detect_degeneracies: " << used_time() - time << "sec \n";

}


// --------------------------------------------------------------------------
// degenerierter incircle test
// --------------------------------------------------------------------------

{

cout << "mache degenerierten incircle Test\n";

real d = sqrt(2);

real a_1=0,b_1= 1,c_1=  d,
     a_2=0,b_2=-1,c_2=  d,
     a_3=1,b_3= 0,c_3= -2*d;

int sign_x=1; int sign_y=1;

real A = a_1*c_2+a_2*c_1; // fr"uher K
real B = b_1*c_2+b_2*c_1; // fr"uher I

real C = 2*c_1*c_2; 

real D = a_1*a_2-b_1*b_2;
real S = a_1*b_2+a_2*b_1;

real RN = sqrt( (a_1*a_1+b_1*b_1)*(a_2*a_2 + b_2*b_2) );

real sqr1 = sqrt(C*(RN + D));
real sqr2 = sqrt(C*(RN - D));

real x_v = A + sign_x*sqr1;
real y_v = B + sign_y*sqr2;
real z_v = RN - (a_1*a_2+b_1*b_2);

real P = a_3*x_v+b_3*y_v+c_3*z_v;

real R = P*P - (a_3*a_3+b_3*b_3)*(x_v*x_v+y_v*y_v);

R.compute_with_precision(60);
R.improve_approximation_to(100);

time = used_time();
assert(R.sign() == 0);
cout << "Zeit fuer degenerierten sign Test: " << used_time() - time << "sec \n";

}


// --------------------------------------------------------------------------
// nicht-degenerierter incircle test
// --------------------------------------------------------------------------

cout << "mache nicht-degenerierten incircle Test\n";

{
real a_1=-12,b_1= -7,c_1=  23,
     a_2=2,b_2=-3,c_2=  9,
     a_3=14,b_3= 10,c_3= 5;

real A = a_1*c_2+a_2*c_1; // fr"uher K
real B = b_1*c_2+b_2*c_1; // fr"uher I

real C = 2*c_1*c_2; 

real D = a_1*a_2-b_1*b_2;

real RN = sqrt( (a_1*a_1+b_1*b_1)*(a_2*a_2 + b_2*b_2) );

real sqr1 = sqrt(C*(RN + D));
real sqr2 = sqrt(C*(RN - D));

real T = abs(a_1)*b_2+abs(a_2)*b_1;
real S = a_1*b_2+a_2*b_1;

int sign_x = T.sign();
int sign_y = sign_x*S.sign();

real x_v = A + sign_x*sqr1;
real y_v = B + sign_y*sqr2;
real z_v = RN - (a_1*a_2+b_1*b_2);

real P = (a_3*x_v+b_3*y_v)+c_3*z_v;
real D3_sq = a_3*a_3+b_3*b_3;
real R_sq = x_v*x_v+y_v*y_v;

real R = P*P - D3_sq*R_sq;

R.compute_with_precision(60);
R.improve_approximation_to(100);

time = used_time();
assert(R.sign() != 0);
cout << "Zeit fuer degenerierten sign Test: " << used_time() - time << "sec \n";
}



// --------------------------------------------------------------------------
// schwieriger nicht-degenerierter incircle test
// --------------------------------------------------------------------------


cout << "mache schwierigen nicht-degenerierten incircle Test\n";

{
real a_1=-12,b_1= -7,c_1=  23,
     a_2=2,b_2=-3,c_2=  9,
     a_3=14,b_3= 10,c_3= 5;

real A = a_1*c_2+a_2*c_1; // fr"uher K
real B = b_1*c_2+b_2*c_1; // fr"uher I

real C = 2*c_1*c_2; 

real D = a_1*a_2-b_1*b_2;

real RN = sqrt( (a_1*a_1+b_1*b_1)*(a_2*a_2 + b_2*b_2) );

real sqr1 = sqrt(C*(RN + D));
real sqr2 = sqrt(C*(RN - D));

real T = abs(a_1)*b_2+abs(a_2)*b_1;
real S = a_1*b_2+a_2*b_1;

int sign_x = T.sign();
int sign_y = sign_x*S.sign();

real x_v = A + sign_x*sqr1;
real y_v = B + sign_y*sqr2;
real z_v = RN - (a_1*a_2+b_1*b_2);

real P = (a_3*x_v+b_3*y_v)+c_3*z_v;
real D3_sq = a_3*a_3+b_3*b_3;
real R_sq = x_v*x_v+y_v*y_v;

real R = P*P - D3_sq*R_sq;

R.compute_with_precision(60);
R.improve_approximation_to(100);

time = used_time();
assert(R.sign() != 0);
cout << "Zeit fuer schwierigen nicht-degenerierten sign Test: " << used_time() - time << "sec \n";
}


// --------------------------------------------------------------------------
//  Oliviers Test mit 2x2 Determinante.
// --------------------------------------------------------------------------

{

cout << "Oliviers Test mit der Determinante\n";
 
real a11 =  (double)  6576833001171515.0;
real a21 =  (double)  4645552482121323.0;
real a12 =  (double) -4901643281669312.0;
real a22 =  (double) -3462280570234346.0;
  
assert((a11*a22-a12*a21).sign() == 1);

a11 = 40000000000000.0; a12 = a11+1; a21 = a11+2; a22 = a11+3;
 
assert((a11*a22-a12*a21).sign() == -1); 

}

// --------------------------------------------------------------------------
//  Stefan Lutters Test mit Wurzeln
// --------------------------------------------------------------------------

{

cout << "Stefan Lutters Test mit der Determinante\n";

time = used_time();
real H = 3*sqrt(real(2))/sqrt(real(20));
real I = (sqrt(real(8)) + sqrt(real(2)))/(2*sqrt(real(5)));
assert((H-I).sign() == 0);
cout << "Zeit fuer Lutter Test: " << used_time() - time << "sec \n";

}


// --------------------------------------------------------------------------
//  Kurts Test mit 3. binomischer Formel.
//  -------------------------------------------------------------------------

{

cout << "Kurts Test mit der 3. binomischen Formel\n";

integer a = 1;
real b,w,z;

cout << "length ";
while (a.length() <= 50)
{
  a = 2*a;
  b = a;
  cout << a.length() << " ";
  w = (sqrt(b+2) - sqrt(b))*(sqrt(b+1) + sqrt(b)) - 1;
  z = (sqrt(b+1) - sqrt(b))*(sqrt(b+1) + sqrt(b)) - 1;
  assert(w.sign() != 0); 
  assert(z.sign() == 0); 
}
cout << "\n";

}

// --------------------------------------------------------------------------
// nochmal Kurts Test mit 3. binomischer Formel, zum Messen der Laufzeit.
// --------------------------------------------------------------------------

{

cout << "Kurts Test mit der 3. binomischen Formel, zur Laufzeitmessung\n";

random_source Source;
int K,L;
int s;
integer Ran;
real b, sqr_b, sqr_b_inc, w;

for (K=1;K<=4000;K*=2)
{
  cout << "k = " << K << "\n";
  b = pow2(integer(K));
  
  if (K <= 31) 
  { 
    Source.set_precision(K);
    Source >> s;
    b = s;
  }
  else
  {
    L = K;
    Ran = 0;
    Source.set_precision(31);
    while(L >= 0)
    {
      Ran = Ran << 32;
      Source >> s;
      Ran += s;
      L = L-32;
    }
    b = Ran;
  }

  sqr_b = sqrt(b);
  sqr_b_inc = sqrt(b+1);
  w = (sqr_b_inc - sqr_b)*(sqr_b_inc + sqr_b);
  time = used_time();
  s = (w-1).sign();
  time = used_time(time);
  cout << " hard sign braucht Zeit " << time << "\n";
  time = used_time();
}

}

// --------------------------------------------------------------------------
//  entarteter Test mit vielen Wurzeln; zum Laufzeitmessen
// ----------------------------------------------------------------------------

{

real y=0, x=0;
for (int ij = 2; ij<= 100; ij++)
{
  x += ij*sqrt(real(ij));
  y = x-x;
  time = used_time();
  assert(y.sign() == 0);
  time = used_time(time);
  cout << ij << " Wurzeln brauchen Zeit " << time << "\n" << flush;
  if (time > 1) break;
}

}

// --------------------------------------------------------------------------
//  beruehmte entartete Tests
// ----------------------------------------------------------------------------


{
  real one = 1, two = 2, three = 3, four = 4, five = 5, six = 6, seven = 7, 
       nine = 9, twelve = 12, twenty = 20, twenty_five = 25, twenty_nine = 29;
  real root_two = sqrt(two); 
  real root_three = sqrt(three);
  real root_five = sqrt(five);
  real root_six = sqrt(six);
  real root_seven = sqrt(seven);
  real third_root_two = root(two,3);
  real third_root_three = root(three,3);
  real third_root_four = root(four,3);
  real third_root_five = root(five,3);
  real third_root_nine = root(nine,3);
  real third_root_twenty = root(twenty,3);
  real third_root_twenty_five = root(twenty_five,3);
  real root_twenty_nine = sqrt(twenty_nine);
  real fifth_root_three = root(three,5);
  real fifth_root_five = root(five,5);
 
  real test;

  time = used_time();
  
  test = sqrt(9+4*root_two)-(1+2*root_two);
  if (sign(test) != 0)
     error_handler(1,"Davenport fails");

  test = sqrt(5+2*root_six) + sqrt(5-2*root_six) - 2*sqrt(three);
  assert(sign(test) == 0);
  // Zippel 

  test = sqrt(16-2*root_twenty_nine+2*sqrt(55-10*root_twenty_nine))
       - (sqrt(22+2*root_five)-sqrt(11+2*root_twenty_nine) + root_five);
  assert(sign(test) == 0);
  // Shanks 

  test = sqrt((112+70*root_two)+(46+34*root_two)*root_five)
       - (5+4*root_two) - (3+root_two)*root_five;
  assert(sign(test) == 0);
  // Borodin

  test = sqrt(3+root_two)*sqrt(6+root_seven)-(3+root_two+root_seven);
  assert(sign(test) != 0);
  // Johannes' Ding ist gar nicht degeneriert.

  time = used_time(time);
  cout << "diverse degenerierte sqrt Tests brauchen Zeit ";
  cout << time << "\n" << flush;

  // ----------------------------------------------------------------------

  time = used_time();

  test = 3*sqrt(third_root_five-third_root_four)
         -third_root_two-third_root_twenty+third_root_twenty_five;
  assert(sign(test) == 0);
  // Ramanujans erster Test wie im Anhang von Johannes' Doktorarbeit
  // deg = 2*3^5 = 486

  test =  root(three,3)*root(seven*third_root_twenty-19,6) 
         -third_root_five + root(two,3);
  assert(sign(test) == 0);
  // Ramanujans zweiter Test wie im Anhang von Johannes' Doktorarbeit,
  // nur divisions-frei
  // deg = 3*3*6*3*3 = 2*3^5 = 486

  test = root(root_five+2,3)-root(root_five-2,3)-1;
  assert(sign(test) == 0);
  // Landaus Test aus SIAM Journal on Computing 19992, p. 96.
  // deg = 2*3*3 = 18

  test = sqrt(3+2*root_three)-one/four*root(twelve,4)*(2+sqrt(twelve));
  assert(sign(test) == 0);
  // Borodin 1 aus Journal of Sumbolic Computation 1985, 196ff.
  // deg = 2*2*4*2 = 2^5 = 32

  test = sqrt(4+3*root_two) - root(two,4)*(1+root_two);
  assert(sign(test) == 0);
  // Borodin 2 aus Journal of Sumbolic Computation 1985, 196ff.
  // deg = 2*2*4 = 16

  time = used_time(time);
  cout << "diverse divisionsfreie Tests mit d <= 6 brauchen Zeit ";
  cout << time << "\n" << flush;

  // ------------------------------------------------------------------------
  
  time = used_time();

  test = root(third_root_two-1,3) -one/third_root_nine +
         root(two,3)/third_root_nine-sqr(third_root_two)/third_root_nine;
  // Ramanujans dritter Test wie in Landaus Artikel,
  // mit Divisionen
  assert(sign(test) == 0);
  // deg = 3*3*3*3 = 3^4 = 81

  test =  root(seven*sqr(third_root_two)*third_root_five-19,6) 
         -third_root_five/third_root_three + third_root_two/third_root_three;
  assert(sign(test) == 0);
  // Ramanujans zweiter Test wie im Anhang von Johannes' Doktorarbeit,
  // mit Divisionen
  // deg = 3*3*6*3 = 2*3^4 = 162

  test = root(two/fifth_root_five
         -fifth_root_three*fifth_root_three*fifth_root_three/fifth_root_five,3)
         -(1+fifth_root_three-sqr(fifth_root_three))/sqr(fifth_root_five);
  assert(sign(test) == 0);
  // Zippels Test (Ramanujan?) 
  // deg = 3^2 * 5*5 = 225

  test = root(three/fifth_root_five
         -fifth_root_three*fifth_root_three*fifth_root_three/fifth_root_five,3)
         -(1+fifth_root_three-sqr(fifth_root_three))/sqr(fifth_root_five);
  assert(sign(test) != 0);
  // Zippel perturbed

  time = used_time(time);
  cout << "Tests mit d <= 6 und Division brauchen Zeit ";
  cout << time << "\n" << flush;


}
  

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

cout << "Alle Tests erfolgreich absolviert\n";

 return 0;

} // main
