//===================================================================
// Schroder's Methods A and B for the calculation of square roots.
// Reference: "On Infinitely Many Algorithms For Solving Equations"
//      by Ernst Schroder, Translated by G. W. Stewart.
// ----------------------------------------------------------------
// For each value of Omega, we calculate A and B then average.
// ----------------------------------------------------------------
// For this family of equations, we consider the equations when the
// parameter Lambda is zero.  Some methods where Lambda is not zero
// will generate identical equations.  Some of these are noted in
// the comments.  For square root, method A and Method B for the
// same Lamda and Omega are closely related to each other.  For
// each equation, if we are finding the square root of c, then our
// estimate k can be used to form the other equation by the
// transformation: y = c/k.  This is also obvious since:
//    c/sqrt( c ) = c^( 1/2 )
// gives us the square root of c also, with an estimate k, then c/k
// gives us a new estimate.  Taking the average is the equivalent
// of a Newton step at the end.
// ----------------------------------------------------------------
// Copyright (C) 1997 by Dann Corbit
//===================================================================
#include <iostream.h>
#include "sqrts.h"
//===================================================================
// Schroder's Method( A, Lamda=0, Omega = 1 ) for x*x - c = 0
// SqrtNewt() produces 2 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtNewt(NUMBER x, NUMBER c)
{
    return (x * x + c) / x * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 1 ) for x*x - c = 0
// c/ans0 is Schroder's Method( B, Lamda=0, Omega = 1 ) for x*x - c = 0
// SqrtAB01() produces 4 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct. Ans0 is also "Newton's Method" for square root.
//===================================================================
NUMBER          SqrtAB01(NUMBER x, NUMBER c)
{
    NUMBER          ans0 = (x * x + c) / x * valOneHalf;
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 2 ) for x*x - c = 0
// Schroder's Method( B, Lamda=-1/2, Omega = 1 ) for x*x - c = 0 and
// Schroder's Method( A, Lamda=1/2, Omega = 1 ) for x*x - c = 0 also
// result in exactly the same equation.
// Interesting, that lower order methods ( Omega=1 ) arrive at the
// same equations as a higher order method ( Omega=2 ). This shows
// that the convergence rate is a minimum limit, and that for a
// given equation, at least one higher order of convergence may be
// achieved. It may or may not be computationally more efficient.
// SqrtAB02() produces 6 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB02(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          ans0 = x * (x2 + val3 * c) / (val3 * x2 + c);
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 3 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega = 3 ) for x*x - c = 0
// and Schroder's Method( B, Lamda=1, Omega = 3 ) for x*x - c= 0
// SqrtAB03() produces 8 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB03(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          ans0 = x - ((x2 - c) * (val3 * x2 + c) * valOneQuarter) /
    (x * (x2 + c));
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 4 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega = 4 ) for x*x - c = 0
// SqrtAB04() produces 10 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB04(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          ans0 = x - val4 * (x2 - c) * x * (x2 + c) /
    (val5 * x2 * (x2 + val2 * c) + c * c);
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 5 ) for x*x - c = 0
// and Schroder's Method( B, Lamda=-1, Omega = 5 ) for x*x-c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega = 5 ) for x*x - c = 0
// and Schroder's Method( A, Lamda=1, Omega = 5 ) for x*x - c = 0
// SqrtAB05() produces 12 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB05(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          ans0 = (x2 + c) * (x2 * (x2 + val14 * c) + c * c) /
    (x * (val3 * x2 + c) * (x2 + val3 * c)) * valOneHalf;
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 6 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega = 6 ) for x*x - c = 0
// SqrtAB06() produces 14 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB06(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          c2 = c * c;
    NUMBER          ans0 = x -
    val2 * (x2 - c) * x * (val3 * (x4 + c2) + val10 * x2 * c) /
    (x4 * (val7 * x2 + val35 * c) + c2 * (val21 * x2 + c));
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 7 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega = 7 ) for x*x - c = 0
// SqrtAB07() produces 16 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB07(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x8 = x4 * x4;
    NUMBER          c2 = c * c;
    NUMBER          c4 = c2 * c2;
    NUMBER          ans0 =
    (val28 * (x2 * c * (x4 + c2)) + (val70 * c2 + x4) * x4 + c2 * c2) /
    (x * (x2 + c) * (x4 + val6 * x2 * c + c2)) * valOneEighth;
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 8 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega = 8 ) for x*x - c = 0
// SqrtAB08() produces 18 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB08(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          c2 = c * c;
    NUMBER          ans0 = x +
    val8 * (c - x2) * x * (c + x2) * (x4 + val6 * x2 * c + c2) /
    ((c + val3 * x2) *
     (val3 * x4 * x2 + val27 * x4 * c + val33 * x2 * c2 + c2 * c));
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega = 9 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega = 9 ) for x*x - c = 0
// SqrtAB09() produces 20 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB09(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x2c = x2 * c;
    NUMBER          c2 = c * c;
    NUMBER          x8 = x4 * x4;
    NUMBER          c4 = c2 * c2;
    NUMBER          ans0 = val2 * x * c * (x4 + val10 * x2c + val5 * c2) * (val5 * x4 + val10 * x2c + c2) /
    ((x2 + c) * (x8 + val44 * x4 * x2 * c + val166 * c2 * x4 + val44 * x2 * c2 * c + c4));
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega =10 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega =10 ) for x*x - c = 0
// SqrtAB010() produces 22 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB10(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x8 = x4 * x4;
    NUMBER          x6 = x4 * x2;
    NUMBER          c2 = c * c;
    NUMBER          c3 = c2 * c;
    NUMBER          c4 = c2 * c2;
    NUMBER          ans0 = x - val2 * (x2 - c) * x * (val5 * (x8 + c4) + val60 * (x6 * c + x2 * c3) + val126 * c2 * x4) /
    (x8 * (val11 * x2 + val165 * c) + val462 * x6 * c2 + val330 * c3 * x4 + c4 * (val55 * x2 + c));
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega =11 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega =11 ) for x*x - c = 0
// SqrtAB011() produces 24 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB11(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x2c = x2 * c;
    NUMBER          c2 = c * c;
    NUMBER          x8 = x4 * x4;
    NUMBER          c4 = c2 * c2;
    NUMBER          ans0 = (x4 + val6 * x2c + c2) * (x8 + val60 * x4 * x2 * c + val134 * c2 * x4 + val60 * x2 * c2 * c + c4) /
    (x * (x2 + c) * (x2 + val3 * c) * (val3 * x2 + c) * (x4 + val14 * x2c + c2)) * valOneQuarter;
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega =12 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega =12 ) for x*x - c = 0
// SqrtAB012() produces 26 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB12(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x8 = x4 * x4;
    NUMBER          x10 = x8 * x2;
    NUMBER          x6 = x4 * x2;
    NUMBER          c2 = c * c;
    NUMBER          c3 = c2 * c;
    NUMBER          c4 = c2 * c2;
    NUMBER          c5 = c4 * c;
    NUMBER          ans0 = x - val4 * (x2 - c) * x * (val3 * x10 + val55 * x8 * c + val198 * x6 * c2 + val198 * c3 * x4 + val55 * x2 * c4 + val3 * c5) /
    (val13 * x8 * x4 + val286 * x10 * c + val1287 * x8 * c2 + val1716 * c3 * x6 + val715 * x4 * c4 + val78 * x2 * c5 + c4 * c2);
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega =13 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega =13 ) for x*x - c = 0
// SqrtAB013() produces 28 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB13(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x8 = x4 * x4;
    NUMBER          c2 = c * c;
    NUMBER          c3 = c2 * c;
    NUMBER          x6 = x4 * x2;
    NUMBER          c4 = c2 * c2;
    NUMBER          x4c = x4 * c;
    NUMBER          x2c2 = x2 * c2;
    NUMBER          ans0 = (x2 + c) * (x8 * x4 + val90 * x8 * x2 * c + val911 * x8 * c2 + val2092 * c3 * x6 + val911 * x4 * c4 + val90 * x2 * c4 * c + c4 * c2) /
    (x * (x6 + val21 * x4c + val35 * x2c2 + val7 * c3) * (val7 * x6 + val35 * x4c + val21 * x2c2 + c3)) * valOneHalf;
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// Ans0 is Schroder's Method( A, Lamda=0, Omega =14 ) for x*x - c = 0
// c / ans0 is Schroder's Method( B, Lamda=0, Omega =14 ) for x*x - c = 0
// SqrtAB014() produces 30 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB14(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x6 = x4 * x2;
    NUMBER          x8 = x4 * x4;
    NUMBER          x14 = x8 * x6;
    NUMBER          c2 = c * c;
    NUMBER          c3 = c2 * c;
    NUMBER          c4 = c2 * c2;
    NUMBER          c7 = c4 * c3;
    NUMBER          x2c6 = x2 * c4 * c2;
    NUMBER          x4c5 = x4 * c4 * c;
    NUMBER          x6c4 = x6 * c4;
    NUMBER          x8c3 = x8 * c3;
    NUMBER          x10c2 = x8 * x2 * c2;
    NUMBER          x12c = x8 * x4 * c;
    NUMBER          ans0 = x * (val105 * x12c + val1365 * x10c2 + val455 * x2c6 + val5005 * x8c3 + val3003 * x4c5 + val6435 * x6c4 + val15 * c7 + x14) /
    (val455 * x12c + val3003 * x10c2 + val105 * x2c6 + val6435 * x8c3 + val1365 * x4c5 + val5005 * x6c4 + c7 + val15 * x14);
    return (ans0 + c / ans0) * valOneHalf;
}

//===================================================================
// c/Ans1 is Schroder's Method( A, Lamda=0, Omega =15 ) for x*x - c = 0
// Ans1 is Schroder's Method( B, Lamda=0, Omega =15 ) for x*x - c = 0
// SqrtAB015() produces 32 times as many correct digits after each
// iteration, if the initial guess x has at least one significant
// digit correct.
//===================================================================
NUMBER          SqrtAB15(NUMBER x, NUMBER c)
{
    NUMBER          x2 = x * x;
    NUMBER          x4 = x2 * x2;
    NUMBER          x6 = x4 * x2;
    NUMBER          x8 = x4 * x4;
    NUMBER          x10 = x8 * x2;
    NUMBER          x12 = x8 * x4;
    NUMBER          x14 = x8 * x6;
    NUMBER          x16 = x8 * x8;
    NUMBER          c2 = c * c;
    NUMBER          c3 = c2 * c;
    NUMBER          c4 = c2 * c2;
    NUMBER          c5 = c4 * c;
    NUMBER          c6 = c4 * c2;
    NUMBER          c7 = c4 * c3;
    NUMBER          c8 = c4 * c4;
    NUMBER          ans1 = val16 * x * c * (val35 * (x12 * c + x2 * c6)
              + val273 * (x10 * c2 + x4 * c5) + val715 * (x8 * c3 + x6 * c4)
                                            + c7 + x14) /
    (val120 * (x14 * c + x2 * c7) + val1820 * (x12 * c2 + x4 * c6)
     + val8008 * (x10 * c3 + x6 * c5) + val12870 * x8 * c4 + x16 + c8);
    return (c / ans1 + ans1) * valOneHalf;
}

#ifdef TEST
#include <iostream.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
void            block(NUMBER x, NUMBER k, int noisy);
void            timeit(NUMBER x, NUMBER k, int which);
int             main(int argc, char **argv)
{
    NUMBER          k;
    NUMBER          x;
#ifdef USE_FLASH
    NUMBER          fudge_factor = "2";
#else
    NUMBER          fudge_factor = 2.;
#endif
    NUMBER          answer;
    long            i;
    long            limit;
    int             noisy = 1;
    int             which;
    char            pszString[256] = "2";
    clock_t         cstart;
    double          delta;
    cout.precision(18);
    if (argc < 2) {
        cout << "Enter a number to find square root of:";
        cin >> pszString;
        cout << pszString;
    } else
        strncpy(pszString, argv[1], sizeof(pszString));
#ifdef USE_FLASH
    k = pszString;
#else
    k = atof(pszString);
#endif
    answer = sqrt(k);
    cout << "\nLib Answer: " << answer;
    for (i = 0; i < 5; i++) {
        fudge_factor = sqrt(fudge_factor);
        x = answer * fudge_factor;      // Guess too high
        if (noisy)
            cout << "\n\nBad Guess High: " << x;
        block(x, k, noisy);
        x = answer / fudge_factor;      // Guess too low
        if (noisy)
            cout << "\n\nBad Guess Low: " << x;
        block(x, k, noisy);
    }
    noisy = 0;
#ifdef USE_FLASH
    limit = 5;
#else
    limit = 500000;
#endif
    for (which = 0; which <= 16; which++) {
        cstart = clock();
        fudge_factor = sqrt(fudge_factor);
        for (i = 0; i < limit; i++) {
            x = answer * fudge_factor;  // Guess too high
            timeit(x, k, which);
            x = answer / fudge_factor;  // Guess too low
            timeit(x, k, which);
            x = answer * fudge_factor;  // Guess too high
            timeit(x, k, which);
            x = answer / fudge_factor;  // Guess too low
            timeit(x, k, which);
            x = answer * fudge_factor;  // Guess too high
            timeit(x, k, which);
            x = answer / fudge_factor;  // Guess too low
            timeit(x, k, which);
        }
        delta = (clock() - cstart) / (double) CLOCKS_PER_SEC;
        if (which == 0)
            cout << "\nNewton-" << which << " " << delta << " seconds.\n" << endl;
        else if (which < 16)
            cout << "\nOmega-" << which << " " << delta << " seconds.\n" << endl;
        else
            cout << "\nEmpty-" << which << " " << delta << " seconds.\n" << endl;
    }
    return 0;
}

void            block(NUMBER x, NUMBER k, int noisy)
{
    NUMBER          answer;
    answer = SqrtAB01(x, k);
    if (noisy)
        cout << "\nAB01 - estimate: " << answer;
    answer = SqrtAB02(x, k);
    if (noisy)
        cout << "\nAB02 - estimate: " << answer;
    answer = SqrtAB03(x, k);
    if (noisy)
        cout << "\nAB03 - estimate: " << answer;
    answer = SqrtAB04(x, k);
    if (noisy)
        cout << "\nAB04 - estimate: " << answer;
    answer = SqrtAB05(x, k);
    if (noisy)
        cout << "\nAB05 - estimate: " << answer;
    answer = SqrtAB06(x, k);
    if (noisy)
        cout << "\nAB06 - estimate: " << answer;
    answer = SqrtAB07(x, k);
    if (noisy)
        cout << "\nAB07 - estimate: " << answer;
    answer = SqrtAB08(x, k);
    if (noisy)
        cout << "\nAB08 - estimate: " << answer;
    answer = SqrtAB09(x, k);
    if (noisy)
        cout << "\nAB09 - estimate: " << answer;
    answer = SqrtAB10(x, k);
    if (noisy)
        cout << "\nAB10 - estimate: " << answer;
    answer = SqrtAB11(x, k);
    if (noisy)
        cout << "\nAB11 - estimate: " << answer;
    answer = SqrtAB12(x, k);
    if (noisy)
        cout << "\nAB12 - estimate: " << answer;
    answer = SqrtAB13(x, k);
    if (noisy)
        cout << "\nAB13 - estimate: " << answer;
    answer = SqrtAB14(x, k);
    if (noisy)
        cout << "\nAB14 - estimate: " << answer;
    answer = SqrtAB15(x, k);
    if (noisy)
        cout << "\nAB15 - estimate: " << answer;
}

void            timeit(NUMBER x, NUMBER k, int which)
{
    NUMBER          answer;
    switch (which) {
    case 0:
        answer = SqrtNewt(x, k);
        break;
    case 1:
        answer = SqrtAB01(x, k);
        break;
    case 2:
        answer = SqrtAB02(x, k);
        break;
    case 3:
        answer = SqrtAB03(x, k);
        break;
    case 4:
        answer = SqrtAB04(x, k);
        break;
    case 5:
        answer = SqrtAB05(x, k);
        break;
    case 6:
        answer = SqrtAB06(x, k);
        break;
    case 7:
        answer = SqrtAB07(x, k);
        break;
    case 8:
        answer = SqrtAB08(x, k);
        break;
    case 9:
        answer = SqrtAB09(x, k);
        break;
    case 10:
        answer = SqrtAB10(x, k);
        break;
    case 11:
        answer = SqrtAB11(x, k);
        break;
    case 12:
        answer = SqrtAB12(x, k);
        break;
    case 13:
        answer = SqrtAB13(x, k);
        break;
    case 14:
        answer = SqrtAB14(x, k);
        break;
    case 15:
        answer = SqrtAB15(x, k);
        break;
    default:
        // Empty timing loop pass
        break;
    }
    return;
}

#endif