#include <math.h>
#include <stdio.h>
#include <stdlib.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));
}
long double     h6(long double x)
{
    return (-120 / (x * x * x * x * x * x));
}

/*******************************/
/* 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(int argc, char **argv)
{
    char           *pszFmtAL1 = "%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             family;
    int             j;
    long double     x;
    long double     bottom = -10;
    long double     top = 10;
    long double     Lamda;
    long double     answer = 0;
    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;

    if (argc < 2)
        family = 0;
    else
        family = atoi(argv[1]);

    k = 4.0;
    for (k = 4.; k <= 400.; k *= 10.) {
        for (fuzz = 1.01; fuzz <= 1.04; fuzz += .01) {

            bottom = -2.5;
            top = 5.;
            switch (family) {
            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;
            default:
                puts("Please enter an option from 0 to 3 on the command line.");  
                puts("\tThe number 0 will analyze square roots from the square function.");
                puts("\tThe number 1 will analyze cubic roots from the cubic function.");
                puts("\tThe number 2 will analyze the exponential from the log function.");
                puts("\tThe number 3 will analyze the log from the exponential function.");
                exit(1);
            }
            for (j = 0; j < 2; j++) {
                if (j == 0) {
                    x = answer * fuzz;
                } else {
                    x = answer / fuzz;
                }
                puts("");
                for (Lamda = bottom; Lamda <= top; Lamda += step) {
                    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);
                    printf(pszFmtAL1, Lamda, errAL3, errBL3);
                }
            }
        }
    }
    return 0;
}