#include <math.h>
#include <stdio.h>
#include <string.h>
#include "schroder.h"
#ifdef DO_ABSOLUTE
#define ABSOLUTE(x) fabsl(x)
#else
#define ABSOLUTE(x) (x)
#endif
/* 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 ) );
}

/*  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;
}
         
/*******************************/
/* 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 = "%Le %Le %Le %Le %Le %Le %Le %Le %Le \n";
   long double ansAL1;
   long double ansAL2;
   long double ansAL3;
   long double ansAL4;
   long double ansBL1;
   long double ansBL2;
   long double ansBL3;
   long double ansBL4;
   long double errAL1;
   long double errAL2;
   long double errAL3;
   long double errAL4;
   long double errBL1;
   long double errBL2;
   long double errBL3;
   long double errBL4;
   int i;
   int j;
   long double x;
   long double bottom = -10;
   long double top = 10;
   long double Lamda;
   long double answer=0;
   long double ffuzz;
   long double minerr;
   long double minlamda=-9999.;
   long double step = 1.0/128.;
   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;


   k = 4.0;
   for ( k = 4.; k <= 400.; k *= 10.)
   {
      for ( fuzz = 1.01; fuzz <= 1.04; fuzz += .01 )
      {
         for ( i = 0; i <= 0; i++ )
         {
            bottom = -2.5;
            top = 5.;
            switch ( i )
            {
               case 0:
                  answer = sqrtl( k );
                  function = f;
                  f1prime  = f1;
                  f2prime  = f2;
                  f3prime  = f3;
                  f4prime  = f3;
                  f5prime  = f3;
                  break;
               case 1:
                  answer = powl( k, 1.0e0L/3.0e0L );
                  function = g;
                  f1prime  = g1;
                  f2prime  = g2;
                  f3prime  = g3;
                  f4prime  = g4;
                  f5prime  = g4;
                  break;
               case 2:
                  answer = expl( k );
                  function = h;
                  f1prime  = h1;
                  f2prime  = h2;
                  f3prime  = h3;
                  f4prime  = h4;
                  f5prime  = h5;
                  break;
               case 3:
                  answer = logl( k );
                  function = o;
                  f1prime  = o1;
                  f2prime  = o1;
                  f3prime  = o1;
                  f4prime  = o1;
                  f5prime  = o1;
                  break;
               case 4:
                  bottom = -50;
                  top = 50;
                  answer = logl( k );
                  function = ldf;
                  f1prime  = ldfPrime;
                  f2prime  = ldf2Prime;
                  f3prime  = ldf3Prime;
                  f4prime  = ldf4Prime;
                  f5prime  = ldf5Prime;
                  break;
            }
            for ( j = 0; j < 2; j++ )
            {
               if ( j == 0 )
               {
                  x = answer * fuzz;
                  ffuzz = fuzz;
               }
               else
               {
                  x = answer / fuzz;
                  ffuzz = 1./fuzz;
               }
               minerr = 1e10;
               for ( Lamda = bottom; Lamda <= top; Lamda += step )
               {
                  ansAL1 = AL1Iteration( x, Lamda, function, f1prime );
                  errAL1= ABSOLUTE( ansAL1-answer );
                  ansBL1 = BL1Iteration( x, Lamda, function, f1prime, f2prime );
                  errBL1= ABSOLUTE( ansBL1-answer );
                  ansAL2 = AL2Iteration( x, Lamda, function, f1prime, f2prime );
                  errAL2= ABSOLUTE( ansAL2-answer );
                  ansBL2 = BL2Iteration( x, Lamda, function, f1prime, f2prime, f3prime );
                  errBL2= ABSOLUTE( ansBL2-answer );
                  ansAL3 = AL3Iteration( x, Lamda, function, f1prime, f2prime, f3prime );
                  errAL3= ABSOLUTE( ansAL3-answer );
                  ansBL3 = BL3Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime );
                  errBL3= ABSOLUTE( ansBL3-answer );
                  ansAL4 = AL4Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime );
                  errAL4= ABSOLUTE( ansAL4-answer );
                  ansBL4 = BL4Iteration( x, Lamda, function, f1prime, f2prime, f3prime, f4prime, f5prime );
                  errBL4= ABSOLUTE( ansBL4-answer );
                  printf( pszFmtAL1, Lamda, errAL1, errBL1, errAL2, errBL2, errAL3, errBL3, errAL4, errBL4 );
               }
            }
         }
      }
   }
   return 0;
}