/*
 * ----------------------------------------------------------------
 * 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 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;
}