#include <stdio.h>
#include <math.h>
#include "Problem6.c"	// Includes the helper functions for the actual Problem
#include "TestCase1.c"	// Includes the helper functions for square root of 2
// A macro to ease printing the roots for the results
#define	PRINT_ROOT(a)	printf("root: %.18Le\n", (long double)a);

// Epsilon as 10^-18
#define eps (long double)0.000000000000000001


// Function Prototypes, so that main could be moved up, works better for
// Commenting all the code.
long double Bisection(long double (*f)(long double), long double a, long double b);
short sign(long double x);
long double Bisection(long double (*f)(long double), long double a, long double b);
long double Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0);
long double Secant(long double (*f)(long double), long double a, long double b);
long double Multi_Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0, int p);


// Main function, calls the necessary functions for the project
// Liberal use of function pointers follows. Quick reference:
// &funtion_name returns a pointer to a functions
// (returntype)(*var)(paramters) is the type for a function pointer saved in var
// (*var)(arguments) calls a function pointer stored in var.
// -----------------------------------------------------------------------------
// Implementations of the three root finding methods follow this function
int main() {
	long double root;
	printf("Numerical Analysis Project 2, Problem 6 Approximator\n");
	printf("Using epsilon = %.11Le\n\n", eps);

	// Calls the test case, from TestCase1.c
	printf("Test Case, 0=x^2-2, guess between 0 and 3, or 1.5 for Newton's method\n");
	root = Bisection(&x_squared_minus_2, 0.0, 3.0);
	PRINT_ROOT(root);
	root = Newton(&x_squared_minus_2, &x_squared_minus_2_prime, 1.5);
	PRINT_ROOT(root);
	root = Secant(&x_squared_minus_2, 0.0, 3.0);
	PRINT_ROOT(root);
	printf("\n");

	// Solves the problem, Problem 6 from the sheet, Electrical Engineering
	printf("Problem 6, unintentional guesses of 0 and 3, or 1.5 for Newton's method\n");
	root = Bisection(&Problem6, 0.0, 3.0);
	PRINT_ROOT(root);
	root = Newton(&Problem6, &Problem6_prime, 1.5);
	PRINT_ROOT(root);
	root = Multi_Newton(&Problem6, &Problem6_prime, 1.5, 2);
	PRINT_ROOT(root);
	root = Secant(&Problem6, 0.0, 3.0);
	PRINT_ROOT(root);
	printf("\n");

	printf("Problem 6, actual guesses of 1*10^10 and 2*10^10, or 1.5*10^10 for Newton\n");
	root = Bisection(&Problem6, 1.0*powl(10.0, 10.0), 2.0*powl(10.0, 10.0));
	PRINT_ROOT(root);
	root = Newton(&Problem6, &Problem6_prime, 1.5*powl(10.0, 10.0));
	PRINT_ROOT(root);
	root = Multi_Newton(&Problem6, &Problem6_prime, 1.5*powl(10.0, 10.0), 2);
	PRINT_ROOT(root);
	root = Secant(&Problem6, 1.0*powl(10.0, 10.0), 2.0*powl(10.0, 10.0));
	PRINT_ROOT(root);
	
	return 0;
}


// Utility function to determine the sign of a long double
// returns 1 for positive or 0 and -1 for < 0
short sign(long double x) {
	if(x < 0.0)
		return -1;
	return 1;
}


// Implementation of the Bisection method
long double Bisection(long double (*f)(long double), long double a, long double b) {
	long double c;	// Center value
	int i = 1;	// Counter value to measure the speed of convergence
	
	// Loop that runs forever, incrementing i on each iteration
	for(;;i++) {
		// Calculate the new center value for this iteration
		c = (a + b) / 2.0;
		// If the difference between b and c is < epsilon we found the root
		if((b-c) < eps) {
			printf("Bisection method in %d iterations\n", i);
			return c;
		}
		// if the root is between b and c
		if( sign(((*f)(c)) * sign((*f)(b))) < 0)
			// set a to be c
			a = c;
		else
			// set b to be c
			b = c;
		// Next iteration
	}
}


// Implementation of Newton's Method
long double Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0) {
	int i = 1;				// Counter value to prevent Newton from running
							// too many iterations
	long double x_1 = 0.0;	// next x value
	long double fprime_of_x_0;	// Derivative of f at the current x value
	
	// Loop that runs forever, incrementing i each time
	for(;;i++) {
		// Calculate the derivative of the current x value
		fprime_of_x_0 = (*fprime)(x_0);
		// If the derivative of f is 0 at the current x value, Newton fails
		if(fabsl(fprime_of_x_0) < eps) {
			printf("WARN: Newton fails, f'(x_0) = 0\n");
			return 0.0;
		}
		
		// Calculate the next x value
		x_1 = x_0 - ((*f)(x_0) / fprime_of_x_0);
		// If the difference between the current and next x values is very small
		// We have found the root.
		if(fabsl(x_1 - x_0) < eps) {
			printf("Newton's method in %d iterations\n", i);
			return x_1;
		}
		// Counter case to ensure Newton's Method doesn't run forever
		if(i > 150) { // Default 50 iterations before failure, bumped up
			printf("WARN: Newton took too long to converge, quitting\n");
			return x_1;
		}
		//Set the current x value equal to the initial for the next iteration.
		x_0 = x_1;
	}
}

// Implementation of Newton's Method for Multiple Roots
long double Multi_Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0, int p) {
	int i = 1;				// Counter value to prevent Newton from running
							// too many iterations
	long double x_1 = 0.0;	// next x value
	long double fprime_of_x_0;	// Derivative of f at the current x value
	
	// Loop that runs forever, incrementing i each time
	for(;;i++) {
		// Calculate the derivative of the current x value
		fprime_of_x_0 = (*fprime)(x_0);
		// If the derivative of f is 0 at the current x value, Newton fails
		if(fabsl(fprime_of_x_0) < eps) {
			printf("WARN: Newton fails, f'(x_0) = 0\n");
			return 0.0;
		}
		
		// Calculate the next x value
		x_1 = x_0 - ((long double)p*(*f)(x_0) / fprime_of_x_0);
		// If the difference between the current and next x values is very small
		// We have found the root.
		if(fabsl(x_1 - x_0) < eps) {
			printf("Newton's method for p=%d in %d iterations\n", p, i);
			return x_1;
		}
		// Counter case to ensure Newton's Method doesn't run forever
		if(i > 150) { // Default 50 iterations before failure, bumped up
			printf("WARN: Newton took too long to converge, quitting\n");
			return x_1;
		}
		//Set the current x value equal to the initial for the next iteration.
		x_0 = x_1;
	}
}


// Implementation of the Secant Method
long double Secant(long double (*f)(long double), long double a, long double b) {
	long double c;	// Center value
	int i = 1;	// Counter value to measure the speed of convergence and
				// prevent secant method from running forever, no guarantee of
				// convergence.
	
	// Loop that runs forever, incrementing i each time
	for(;;i++) {
		// Calculate the center value by the secant method
		c = b - ((long double)(b - a)*(*f)(b)) / ((*f)(b) - (*f)(a));
		// If the difference between the center and the upper bound is small
		// We have found the root
		if(fabsl(c - b) < eps) {
			printf("Secant method in %d iterations\n", i);
			return c;
		}
		// Counter case to ensure the Secant Method doesn't run forever
		if(i > 50) { // Default 50 iterations before failure, bumped up
			printf("WARN: Secant took too long to converge, quitting\n");
			return c;
		}
		// Set the new bounds for the next iteration
		a = b;
		b = c;
	}
}