#include "schroder.h"
/********************************************************************/
/*  Schroder's methods for root calculation, implemented in C.      */
/*  This is the 'generic' version for any function with derivatives */
/*  ----------------------------------------------------------------*/
/*  Reference: "On Infinitely Many Algorithms For Solving Equations"*/
/*             by Ernst Schroder, Translated by G. W. Stewart.      */
/*  ----------------------------------------------------------------*/
/*  Copyright (C) 1996 by Dann Corbit                               */
/* =================================================================*/
/*  Inspired by the following reference:                            */
/*  Local convergence rates can be made arbitrarily large.          */
/*  Ernst Schroder showed, over 100 years ago, how a family         */
/*  of explicit single-point formulas could be designed to          */
/*  have arbitrarily high orders of (local) convergence.            */
/*  You can obtain a translation of Schroder's paper at:            */
/*                                                                  */
/*     ftp://ftp.cs.umd.edu/pub/papers/papers/2990/2990.ps.Z        */
/*                                                                  */
/* =================================================================*/
/*  A posthumous thanks to Mr. Schroder, & thanks to Pete Stewart.  */
/********************************************************************/
/***********************************************************/
/* Schroder's Method, Series A, Lamda variable, Omega = 1. */
/***********************************************************/
long double AL1Iteration(
   long double x,                          /* The input value x   */
   long double L,                          /* The parameter Lamda */
   long double ( *f )( long double ),      /* The function f(x)   */
   long double ( *f1 )( long double )      /* The 1st derivative  */
)
{
   long double Numerator;
   long double Denominator;
   long double f_0 = f( x );
   long double f_1 = f1( x );
   Numerator   = x * f_0;
   Denominator = -f_1 * x + L * f_0;
   return x + Numerator/Denominator;
}

/***********************************************************/
/* Schroder's Method, Series A, Lamda variable, Omega = 2. */
/***********************************************************/
long double AL2Iteration(
   long double x,                       /* The input value x   */
   long double L,                       /* The parameter Lamda */
   long double ( *f )( long double ),   /* The function f(x)   */
   long double ( *f1 )( long double ),  /* The 1st derivative  */
   long double ( *f2 )( long double )   /* The 2nd derivative  */
)
{
   long double Numerator;
   long double Denominator;
   long double f_0 = f( x );
   long double f_1 = f1( x );
   long double f_2 = f2( x );
   Numerator   =  2 * f_0 * x * ( f_1 * x - L * f_0 );
   Denominator = -2 * ( f_1 * f_1 ) * x * x + f_0 * f_2 * x * x
                 +2 * L * f_0 * f_1 * x - L * L * f_0 * f_0 + L * f_0 * f_0;
   return x + Numerator/Denominator;
}

/***********************************************************/
/* Schroder's Method, Series A, Lamda variable, Omega = 3. */
/***********************************************************/
long double AL3Iteration(
   long double x,                       /* The input value x   */
   long double L,                       /* The parameter Lamda */
   long double ( *f )( long double ),   /* The function f(x)   */
   long double ( *f1 )( long double ),  /* The 1st derivative  */
   long double ( *f2 )( long double ),  /* The 2nd derivative  */
   long double ( *f3 )( long double )   /* The 3rd derivative  */
)
{
   long double Numerator;
   long double Denominator;
   long double f_0 = f( x );
   long double f_1 = f1( x );
   long double f_2 = f2( x );
   long double f_3 = f3( x );
   Numerator   =  3 * f_0 * x * (
                    - 2 * f_1 * f_1 * x * x
                    + f_0 * f_2 * x * x
                    + 2 * L * f_0 * f_1 * x
                    - L * L * f_0 * f_0
                    + L * f_0 * f_0 );
   Denominator =  6 * f_1 * f_1 * f_1 * x * x * x
                 -6 * f_0 * f_2 * f_1 * x * x * x
                 +    f_0 * f_0 * f_3 * x * x * x
                 -6 * L * f_0 * f_1 * f_1 * x * x
                 +3 * L * f_0 * f_0 * f_2 * x * x
                 -3 * L * f_0 * f_0 * f_1 * x
                 -    L * L * L * f_0 * f_0 * f_0
                 +3 * L * L * f_0 * f_0 * f_0
                 -2 * L * f_0 * f_0 * f_0;

   return x + Numerator/Denominator;
}

/***********************************************************/
/* Schroder's Method, Series A, Lamda variable, Omega = 4. */
/***********************************************************/
long double AL4Iteration(
   long double x,                       /* The input value x   */
   long double L,                       /* The parameter Lamda */
   long double ( *f )( long double ),   /* The function f(x)   */
   long double ( *f1 )( long double ),  /* The 1st derivative  */
   long double ( *f2 )( long double ),  /* The 2nd derivative  */
   long double ( *f3 )( long double ),  /* The 3rd derivative  */
   long double ( *f4 )( long double )   /* The 4th derivative  */
)
{
   long double t1 = f(x);
   long double t3 = f1(x);
   long double t4 = t3 * t3;
   long double t6 = x * x;
   long double t7 = t6 * x;
   long double t8 = t4 * t3 * t7;
   long double t9 = f2(x);
   long double t10 = t1 * t9;
   long double t13 = t1 * t1;
   long double t14 = f3(x);
   long double t15 = t13 * t14;
   long double t17 = L * t1;
   long double t18 = t4 * t6;
   long double t20 = L * t13;
   long double t21 = t9 * t6;
   long double t23 = L * L;
   long double t24 = t23 * t13;
   long double t25 = t3 * x;
   long double t28 = t23 * L;
   long double t29 = t13 * t1;
   long double t30 = t28 * t29;
   long double t31 = t23 * t29;
   long double t32 = L * t29;
   long double t34 = t4 * t4;
   long double t35 = t6 * t6;
   long double t41 = t9 * t9;
   long double t60 = t23 * t23;
   long double t61 = t13 * t13;
   long double t66 = 24.0 * t34 * t35 - 36.0 * t10 * t4 * t35 + 8.0 * t15 * t3 * t35 + 6.0 * t13 * t41 * t35 - t29 *
      f4(x) * t35 - 24.0 * t17 * t8 + 24.0 * t20 * t9 * t3 * t7 - 4.0 * t32 *
      t14 * t7 + 12.0 * t24 * t18 - 6.0 * t31 * t21 - 12.0 * t20 * t18 + 6.0 * t32 * t21 - 4.0 * t30 * t25 + 12.0 * t31 *
      t25 - 8.0 * t32 * t25 + t60 * t61 - 6.0 * t28 * t61 + 11.0 * t23 * t61 - 6.0 * L * t61;
   long double t70 = x - 4.0 * t1 * x * (6.0 * t8 - 6.0 * t10 * t3 * t7 + t15 * t7 - 6.0 * t17 * t18 + 3.0 * t20 * t21 + 3.0
      * t24 * t25 - 3.0 * t20 * t25 - t30 + 3.0 * t31 - 2.0 * t32) / t66;
   return t70;
}

/***********************************************************/
/* Schroder's Method, Series B, Lamda variable, Omega = 1. */
/***********************************************************/
long double BL1Iteration(
   long double x,                       /* The input value x   */
   long double L,                       /* The parameter Lamda */
   long double ( *f )( long double ),   /* The function f(x)   */
   long double ( *f1 )( long double ),  /* The 1st derivative  */
   long double ( *f2 )( long double )   /* The 2nd derivative  */
)
{
   long double Numerator;
   long double Denominator;
   long double f_0 = f( x );
   long double f_1 = f1( x );
   long double f_2 = f2( x );
   Numerator   = x * f_1 * f_0;
   Denominator = - f_1 * f_1 * x
                 + f_0 * f_2 * x
                 + L * f_0 * f_1;
   return x + Numerator/Denominator;
}

/***********************************************************/
/* Schroder's Method, Series B, Lamda variable, Omega = 2. */
/***********************************************************/
long double BL2Iteration(
   long double x,                       /* The input value x   */
   long double L,                       /* The parameter Lamda */
   long double ( *f )( long double ),   /* The function f(x)   */
   long double ( *f1 )( long double ),  /* The 1st derivative  */
   long double ( *f2 )( long double ),  /* The 2nd derivative  */
   long double ( *f3 )( long double )   /* The 3rd derivative  */
)
{
   long double Numerator;
   long double Denominator;
   long double f_0 = f( x );
   long double f_1 = f1( x );
   long double f_2 = f2( x );
   long double f_3 = f3( x );
   Numerator   = 2 * f_0 * x * (
                  - f_1 * f_1 * x
                  + f_0 * f_2 * x
                  + L * f_0 * f_1 );
   Denominator =  2 * f_1 * f_1 * f_1 * x * x
                 -3 * f_0 * f_2 * f_1 * x * x
                 +    f_0 * f_0 * f_3 * x * x
                 -2 * L * f_0 * f_1 * f_1 * x
                 +2 * L * f_0 * f_0 * f_2 * x
                 +    L * L * f_0 * f_0 * f_1
                 -    L * f_0 * f_0 * f_1;
   return x + Numerator/Denominator;
}

/***********************************************************/
/* Schroder's Method, Series B, Lamda variable, Omega = 3. */
/***********************************************************/
long double BL3Iteration(
   long double x,                       /* The input value x   */
   long double L,                       /* The parameter Lamda */
   long double ( *f )( long double ),   /* The function f(x)   */
   long double ( *f1 )( long double ),  /* The 1st derivative  */
   long double ( *f2 )( long double ),  /* The 2nd derivative  */
   long double ( *f3 )( long double ),  /* The 3rd derivative  */
   long double ( *f4 )( long double )   /* The 4th derivative  */
)
{
   long double Numerator;
   long double Denominator;
   long double f_0 = f( x );
   long double f_1 = f1( x );
   long double f_2 = f2( x );
   long double f_3 = f3( x );
   long double f_4 = f4( x );
   Numerator   = 3 * f_0 * x * (
                    -2 * f_1 * f_1 * f_1 * x * x
                    +3 * f_0 * f_2 * f_1 * x * x
                    -    f_0 * f_0 * f_3 * x * x
                    +2 * L * f_0 * f_1 * f_1 * x
                    -2 * L * f_0 * f_0 * f_2 * x
                    -    L * L * f_0 * f_0 * f_1
                    +    L * f_0 * f_0 * f_1 );
   Denominator =  6  * f_1 * f_1 * f_1 * f_1 * x * x * x
                 -12 * f_0 * f_2 * f_1 * f_1 * x * x * x
                 +4  * f_0 * f_0 * f_3 * f_1 * x * x * x
                 +3  * f_0 * f_0 * f_2 * f_2 * x * x * x
                 -     f_0 * f_0 * f_0 * f_4 * x * x * x
                 -6  * L * f_0 * f_1 * f_1 * f_1 * x * x
                 +9  * L * f_0 * f_0 * f_2 * f_1 * x * x
                 -3  * L * f_0 * f_0 * f_0 * f_3 * x * x
                 +3  * L * L * f_0 * f_0 * f_1 * f_1 * x
                 -3  * L * L * f_0 * f_0 * f_1 * f_1 * x
                 +3  * L * f_0 * f_0 * f_0 * f_2 * x
                 -     L  * L * L * f_0 * f_0 * f_0 * f_1
                 +3  * L * L * f_0 * f_0 * f_0 * f_1
                 -2  * L * f_0 * f_0 * f_0 * f_1;

   return x + Numerator/Denominator;
}

/***********************************************************/
/* Schroder's Method, Series B, Lamda variable, Omega = 4. */
/***********************************************************/
long double BL4Iteration(
   long double x,                       /* The input value x   */
   long double L,                       /* The parameter Lamda */
   long double ( *f )( long double ),   /* The function f(x)   */
   long double ( *f1 )( long double ),  /* The 1st derivative  */
   long double ( *f2 )( long double ),  /* The 2nd derivative  */
   long double ( *f3 )( long double ),  /* The 3rd derivative  */
   long double ( *f4 )( long double ),  /* The 4th derivative  */
   long double ( *f5 )( long double )   /* The 5th derivative  */
)
{
   long double t1 = f(x);
   long double t3 = f1(x);
   long double t4 = t3 * t3;
   long double t5 = t4 * t4;
   long double t6 = x * x;
   long double t7 = t6 * x;
   long double t8 = t5 * t7;
   long double t9 = f2(x);
   long double t10 = t1 * t9;
   long double t13 = t1 * t1;
   long double t14 = f3(x);
   long double t15 = t13 * t14;
   long double t18 = t9 * t9;
   long double t21 = t13 * t1;
   long double t22 = f4(x);
   long double t23 = t21 * t22;
   long double t25 = L * t1;
   long double t26 = t4 * t3;
   long double t27 = t26 * t6;
   long double t29 = L * t13;
   long double t31 = t9 * t3 * t6;
   long double t33 = L * t21;
   long double t34 = t14 * t6;
   long double t36 = L * L;
   long double t37 = t36 * t13;
   long double t38 = t4 * x;
   long double t40 = t36 * t21;
   long double t41 = t9 * x;
   long double t45 = t36 * L;
   long double t46 = t45 * t21;
   long double t50 = 6.0 * t8 - 12.0 * t10 * t4 * t7 + 4.0 * t15 * t3 * t7 + 3.0 * t13 * t18 * t7 - t23 * t7 - 6.0 * t25 *
      t27 + 9.0 * t29 * t31 - 3.0 * t33 * t34 + 3.0 * t37 * t38 - 3.0 * t40 * t41 - 3.0 * t29 * t38 + 3.0 * t33 * t41 - t46
      * t3 + 3.0 * t40 * t3 - 2.0 * t33 * t3;
   long double t51 = t13 * t13;
   long double t52 = L * t51;
   long double t54 = t36 * t51;
   long double t56 = t36 * t36;
   long double t59 = t45 * t51;
   long double t63 = t6 * t6;
   long double t78 = 6.0 * t52 * t3 - 11.0 * t54 * t3 - t56 * t51 * t3 + 6.0 * t59 * t3 - t51 *
      f5(x) * t63 + 12.0 * t29 * t27 - 6.0 * t54 * t34 - 18.0 * t33 * t31 - 4.0 * t59
      * t41 + 4.0 * t46 * t38 - 24.0 * t5 * t3 * t63 - 12.0 * t37 * t27 - 4.0 * t52 * t22 * t7 + 16.0 * t33 * t14 * t3 * t7;
   long double t103 = 24.0 * t25 * t8 + 12.0 * t33 * t18 * t7 - 8.0 * t52 * t41 + 8.0 * t33 * t38 + 12.0 * t54 * t41
      - 12.0 * t40 * t38 + 10.0 * t21 * t9 * t14 * t63 - 30.0 * t3 * t13 * t18 * t63 - 48.0 * t29 * t9 * t4 * t7 - 20.0 *
      t15 * t4 * t63 + 60.0 * t10 * t26 * t63 + 5.0 * t23 * t3 * t63 + 6.0 * t52 * t34 + 18.0 * t40 * t31;
   long double t108 = x + 4.0 * t1 * x * t50 / (t78 + t103);
   return t108;
}