// Arup Guha
// 2/2/2012
// Example of Newton's Algorithm, adjusted for points where the derivative is 0 at the
// point of interest. The function used here is the one from problem 2.4 #2b.

#include <math.h>

double f(double x);
double fprime(double x);
double adjustedNewton(double p0, int maxiter);
double fdoubleprime(double x);

int main() {

    printf("root is %lf\n", adjustedNewton(-3,10000));
    return 0;
}

// The specific function to test the algorithm.
double f(double x) {
    return pow(x,6) + 6*pow(x,5)+9*pow(x,4)-2*x*x*x - 6*x*x + 1;
}

double fprime(double x) {
    return 6*pow(x,5)+30*x*x*x*x+36*x*x*x - 6*x*x - 12*x;
}

double fdoubleprime(double x) {
    return 30*x*x*x*x + 120*x*x*x + 108*x*x - 12*x - 12;
}

// Uses the initial guess init and runs the fixed point
// algorithm for a maximum of maxiter iterations.
double adjustedNewton(double p0, int maxiter) {

    int i=1;

    while (i <= maxiter) {

        // Avoid divide by zero.
        if (fprime(p0) == 0)
            p0 += .01;

        double der_0 = f(p0);
        double der_1 = fprime(p0);
        double der_2 = fdoubleprime(p0);

        double p1 = p0 - (der_0*der_1)/(der_1*der_1-der_0*der_2);

        printf("Iter %d: %.9lf\n", i, p1);

        if (fabs(p1 - p0) < 0.00001)
            return p1;

        i++;
        p0 = p1;
    }

    // Algorithm Failed.
    return -1;

}