#include <math.h>
#include <stdio.h>
#include <string.h>
#include "schroder.h"

/* Some function prototypes for testing: */
long double k = 5;
/*******************************/
/* function family for testing */
/* the square root of x:       */
/*******************************/
long double f( long double x )
{
   return ( x*x - k );
}
long double f1( long double x )
{
   return ( x + x );
}
long double f2( long double x )
{
   return ( 2. );
}
long double f3( long double x )
{
   return ( 0. );
}
/*******************************/
/* function family for testing */
/* the cube root of x:         */
/*******************************/
long double g( long double x )
{
   return ( x*x*x - k );
}
long double g1( long double x )
{
   return ( 3 * x * x );
}
long double g2( long double x )
{
   return ( 6 * x );
}
long double g3( long double x )
{
   return ( 6. );
}
long double g4( long double x )
{
   return ( 0. );
}
/*******************************/
/* function family for testing */
/* the exponential of x:       */
/*******************************/
long double h( long double x )
{
   return ( logl( x ) - k );
}
long double h1( long double x )
{
   return ( 1/x  );
}
long double h2( long double x )
{
   return ( -1/( x*x ) );
}
long double h3( long double x )
{
   return ( 2/( x*x*x ) );
}
long double h4( long double x )
{
   return ( -6/( x*x*x*x ) );
}
long double h5( long double x )
{
   return ( 24/( x*x*x*x*x ) );
}
long double h6( long double x )
{
   return ( -120/( x*x*x*x*x*x ) );
}
/*  Function ldf of x is x^x */
long double ldf( long double x )
{
   return powl( x, x ) - k;
}
/*  Function ldfPrime of x is 1st derivative of x^x */
long double ldfPrime( long double x )
{
   return powl( x, x ) * ( logl( x ) + 1 );
}
/*  Function ldf2Prime of x is 2nd derivative of x^x */
long double ldf2Prime( long double x )
{
   long double xx   = powl( x, x );
   long double xxm1 = powl( x, x - 1 );
   long double lgx  = logl( x );
   long double lgx2 = lgx * lgx;
   return xx * lgx2 + 2 * xx * lgx + xx + xxm1;
}
/*  Function ldf3Prime of x is 3rd derivative of x^x */
long double ldf3Prime( long double x )
{
   long double xx   = powl( x, x );
   long double xxm1 = powl( x, x - 1 );
   long double xxm2 = powl( x, x - 2 );
   long double lgx  = logl( x );
   long double lgx2 = lgx * lgx;
   return xx * lgx2 * lgx + 3 * xx * lgx2 + 3 * xx * lgx + xx + 3 * xxm1 * lgx + 3 * xxm1 - xxm2;
}
long double ldf4Prime( long double x )
{
   long double xx   = powl( x, x );
   long double xxm1 = powl( x, x - 1 );
   long double xxm2 = powl( x, x - 2 );
   long double xxm3 = powl( x, x - 3 );
   long double lgx  = logl( x );
   long double lgx2 = lgx * lgx;
   long double lgx4 = lgx2 * lgx2;
   return xx * lgx4 + 4 * xx * lgx2 * lgx + 6 * xx * lgx2 + 4 * xx * lgx + xx + 6 * xxm1 * lgx2 + 12 * xxm1 * lgx + 6 * xxm1 - xxm2 - 4 * xxm2 * lgx + 2 * xxm3;
}

long double ldf5Prime( long double x )
{
   long double xx = powl( x, x );
   long double lgx2 = logl( x );
   long double lgx4 = lgx2 * lgx2;
   long double lgx6 = lgx4 * lgx2;
   long double lgx8 = lgx6 * lgx2;
   long double xxm1 = powl( x, x - 1 );
   long double xxm2 = powl( x, x - 2 );
   long double xxm3 = powl( x, x - 3 );
   long double xxm4 = powl( x, x - 4 );
   return xx * lgx8 * lgx2 + 5 * xx * lgx8 + 10 * xx * lgx6 + 10 * xx * lgx4 + 5 * xx * lgx2 + xx + 10 * xxm1 * lgx6 + 30 * xxm1 * lgx4
      + 30 * xxm1 * lgx2 + 10 * xxm1 - 5 * xxm2 * lgx2 + 5 * xxm2 - 10 * xxm2 * lgx4 + 10 * xxm3 * lgx2 - 6.0 * xxm4;
}

long double ldf6Prime( long double x )
{
   long double xx   = powl( x, x );
   long double xxm1 = powl( x, x - 1 );
   long double xxm2 = powl( x, x - 2 );
   long double xxm3 = powl( x, x - 3 );
   long double xxm4 = powl( x, x - 4 );
   long double xxm5 = powl( x, x - 5 );
   long double lgx  = logl( x );
   long double lgx2 = lgx * lgx;
   long double lgx4 = lgx2 * lgx2;
   long double lgx3 = lgx2 * lgx;
   long double temp = 60 * xxm1 * lgx3 + 90 * xxm1 * lgx2 + 60 * xxm1 * lgx + 30 * xxm2 * lgx - 15 * xxm2 * lgx2
      + xx * lgx4 * lgx2 + 15 * xxm1 * lgx4 - 20 * xxm2 * lgx3 + 30 * xxm3 * lgx2 - 36 * xxm4 * lgx + 24 * xxm5;
   return 4 * xxm4 - 15 * xxm3 + 25 * xxm2 + 15 * xxm1 + xx + 6 * xx * lgx4 * lgx + 15 * xx * lgx4
      + 20 * xx * lgx3 + 15 * xx * lgx2 + 6 * xx * lgx + temp;
}
/*******************************/
/* function family for testing */
/* the natural log of x:       */
/*******************************/
long double o( long double x )
{
   return ( expl( x ) - k );
}
long double o1( long double x )
{
   return ( expl( x ) );
}

extern int __cdecl main();

/* Test driver routine: */
int main()
{
   char *pszFmtAL1 = "AL1 %.10Lf error: %.10Lf\n";
   char *pszFmtAL2 = "AL2 %.10Lf error: %.10Lf\n";
   char *pszFmtAL3 = "AL3 %.10Lf error: %.10Lf\n";
   char *pszFmtAL4 = "AL4 %.10Lf error: %.10Lf\n";
   char *pszFmtAL5 = "AL5 %.10Lf error: %.10Lf\n";
   char *pszFmtAL6 = "AL6 %.10Lf error: %.10Lf\n";

   char *pszFmtBL1 = "BL1 %.10Lf error: %.10Lf\n";
   char *pszFmtBL2 = "BL2 %.10Lf error: %.10Lf\n";
   char *pszFmtBL3 = "BL3 %.10Lf error: %.10Lf\n";
   char *pszFmtBL4 = "BL4 %.10Lf error: %.10Lf\n";
   char *pszFmtBL5 = "BL5 %.10Lf error: %.10Lf\n";
   char *pszFmtBL6 = "BL6 %.10Lf error: %.10Lf\n";
   int i;
   int j;
   long double x;
   long double xi;
   long double bottom = -10;
   long double top = 10;
   long double Lamda;
   long double answer=0;
   long double error;
   long double minerr;
   long double minlamda=-9999.;
   long double step = 1.0/255.;
   long double fuzz;
   long double ( __cdecl *function ) ( long double ) = f;
   long double ( __cdecl *f1prime ) ( long double ) = f1;
   long double ( __cdecl *f2prime ) ( long double ) = f2;
   long double ( __cdecl *f3prime ) ( long double ) = f3;
   long double ( __cdecl *f4prime ) ( long double ) = f3;
   long double ( __cdecl *f5prime ) ( long double ) = f3;
   long double ( __cdecl *f6prime ) ( long double ) = f3;


   for ( k = 2.0; k <= 5; k+= 0.5 )
   {
      printf( "Roots for %30.20LeL\n", k );
      for ( fuzz = 1.1; fuzz >= 1.00001; fuzz = sqrtl(fuzz) )
      {
         printf( "Fuzz factor %30.20LeL\n", fuzz );
         for ( i = 0; i < 5; i++ )
         {
            bottom = -3;
            top = 3;
            switch ( i )
            {
               case 0:
                  puts( "\nSquare root tests:" );
                  answer = sqrtl( k );
                  function = f;
                  f1prime  = f1;
                  f2prime  = f2;
                  f3prime  = f3;
                  f4prime  = f3;
                  f5prime  = f3;
                  f6prime  = f3;
                  break;
               case 1:
                  puts( "\nCube root tests:" );
                  answer = powl( k, 1.0e0L/3.0e0L );
                  function = g;
                  f1prime  = g1;
                  f2prime  = g2;
                  f3prime  = g3;
                  f4prime  = g4;
                  f5prime  = g4;
                  f6prime  = g4;
                  break;
               case 2:
                  puts( "\nExponential tests:" );
                  answer = expl( k );
                  function = h;
                  f1prime  = h1;
                  f2prime  = h2;
                  f3prime  = h3;
                  f4prime  = h4;
                  f5prime  = h5;
                  f6prime  = h6;
                  break;
               case 3:
                  puts( "\nNatural log tests:" );
                  answer = logl( k );
                  function = o;
                  f1prime  = o1;
                  f2prime  = o1;
                  f3prime  = o1;
                  f4prime  = o1;
                  f5prime  = o1;
                  f6prime  = o1;
                  break;
               case 4:
                  bottom = -50;
                  top = 50;
                  puts( "\nPerfect root tests:" );
                  answer = logl( k );
                  function = ldf;
                  f1prime  = ldfPrime;
                  f2prime  = ldf2Prime;
                  f3prime  = ldf3Prime;
                  f4prime  = ldf4Prime;
                  f5prime  = ldf5Prime;
                  f6prime  = ldf6Prime;
                  break;
            }
            for ( j = 0; j < 2; j++ )
            {
               if ( j == 0 )
               {
                  puts( "overestimate" );
                  x = answer * fuzz;
               }
               else
               {
                  puts( "underestimate" );
                  x = answer / fuzz;
               }
               printf( "Initial guess: %.10Lf\n", x );
               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = AL1Iteration( x, Lamda, function, f1prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtAL1, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = BL1Iteration( x, Lamda, function, f1prime, f2prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtBL1, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = AL2Iteration( x, Lamda, function, f1prime, f2prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtAL2, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = BL2Iteration( x, Lamda, function, f1prime, f2prime, f3prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtBL2, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = AL3Iteration( x, Lamda, function, f1prime, f2prime, f3prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtAL3, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = BL3Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtBL3, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = AL4Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtAL4, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = BL4Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime, f5prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
/*
               printf( pszFmtBL4, minlamda, minerr );

               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = AL5Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime, f5prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }

               printf( pszFmtAL5, minlamda, minerr );
*/
               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  xi = BL5Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime, f5prime, f6prime );
                  error = fabsl( xi-answer );
                  if ( error < minerr )
                  {
                     minlamda = Lamda;
                     minerr = error;
                  }
               }
               printf( pszFmtBL5, minlamda, minerr );
            }
         }
      }
   }
   return 0;
}