#include <math.h>

//===================================================================
// 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.
//===================================================================
long double     SqrtAB02(long double x, long double c)
{
    long double     x2 = x * x;
    long double     ans0 = x * (x2 + 3. * c) / (3. * x2 + c);
    return (ans0 + c / ans0) * 0.5;
}
long double     SqrtGuess(long double d)
{
    long double     z;
    int             e;

    static const long double c0 = 0.4173075996388649989088799756398314815405054464951;
    static const long double c1 = 0.5901620670906445829966311803710009743270304792930;
    static const long double c2 = 0.8346151992777299978177599512796629630810108929897;
    static const long double c3 = 0.2950810335453222914983155901855004871635152396465;

    if (d != 0) {
        z = frexpl(d, &e);
        if (!(e & 1)) {         // exponent is even
            d = c0 + c1 * z;
        } else {                // exponent is odd
            if (e < 0)
                d = c3 + c0 * z;// exponent is odd & negative
            else
                d = c1 + c2 * z;// exponent is odd & positive
        }
        d = ldexpl(d, e / 2);
    }
    return d;
}
#ifdef TEST
#include <stdio.h>
int             main()
{
    char            string[255];
    long double     c;
    long double     x;
    puts("enter a long double to find the square root of:");
    if (fgets(string, sizeof(string), stdin) != NULL) {
        sscanf(string, "%Lf", &c);

        x = SqrtGuess(c);

        printf("Primitive guess is:%f\n", x);
        if (x != 0) {
            printf("square root of %Lg is approximately %.18Lg\n", c, SqrtAB02(x, c));
        }
    }
    return 0;
}
#endif