/*-----------------------------------------------------------------------------
Filename   : root.c
Name       : Tristan Miller
SID#       : 123456789
Description: Finds the root of a function using various numerical methods
-----------------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "function.h"

double bisection(double a, double b, double tol, int max_iter, double (*f)());
double newton_raphson(double p0, double tol, int max_iter,
	      double (*f)(), double(*fp)());
double secant(double p0, double p1, double tol, int max_iter, double (*f)());

/*-----------------------------------------------------------------------------
                                      main()
-----------------------------------------------------------------------------*/
int main(void) {

  double r, tol, p[2];
  int choice, max_iter, i;

  printf("This program finds the root of f(x).\n");

  do {
    printf("Enter the maximum number of iterations: ");
    if (scanf("%d", &max_iter) == EOF) {
      printf("Error -- unexpected end of input\n");
      exit(EXIT_FAILURE);
    }
    if (max_iter <= 0)
      printf("You must iterate at least once!\n");
  } while (max_iter <= 0);

  do {
  printf("Enter the tolerance: ");
    if (scanf("%lf", &tol) == EOF) {
      printf("Error -- unexpected end of input\n");
      exit(EXIT_FAILURE);
    }
    if (tol < 0)
      printf("Tolerance must be non-negative.\n");
  } while (tol < 0);

  for (i = 0; i < 2; i++) {
    printf("Enter the initial guess p%d: ", i);
    if (scanf("%lf", &(p[i])) == EOF) {
      printf("Error -- unexpected end of input\n");
      exit(EXIT_FAILURE);
    }
  }

  while (1) {
    printf("\nMENU\n"
	   "====\n"
	   "1. Bisection method\n"
	   "2. Newton--Raphson method\n"
	   "3. Secant method\n"
	   "4. Quit\n"
	   "Your selection: ");

    choice = 4;
    scanf("%d", &choice);

    switch(choice) {

    case 1:
      printf("Bisection method:\n");
      r = bisection(p[0], p[1], tol, max_iter, f);
      break;

    case 2:
      printf("Newton--Raphson method:\n");
      r = newton_raphson(p[0], tol, max_iter, f, fprime);
      break;

    case 3:      
      printf("Secant method:\n");
      r = secant(p[0], p[1], tol, max_iter, f);
      break;

    case 4:
      return EXIT_SUCCESS;

    default:
      printf("That is not a valid option.\n");
      continue;
    }

    if (r != HUGE_VAL)
      printf("The root is %f.\n", r);
    else
      printf("Root not found.");
  }

}


/*-----------------------------------------------------------------------------
Name   : bisection()
Purpose: finds a root using bisection method
Returns: the root, or HUGE_VAL if no root was found
Args   : a, b     - initial guesses
         tol      - tolerance
         max_iter - maximum number of iterations
         f        - the function
Pre    : f(a)f(b) < 0
-----------------------------------------------------------------------------*/
double bisection(double a, double b, double tol, int max_iter, double (*f)()) {
  double c, increment, result;
  int i = 0;

  printf("Iteration\tApproximation\n"
	 "---------\t-------------\n");

  while (i++ < max_iter) {
    increment = (b - a) / 2.0;
    c = a + increment;
    result = f(c);

    printf("%9d\t%.12f\n", i, c);

    if (increment < tol || result == 0.0)
      return c;

    if (f(a) * result > 0)
      a = c;
    else
      b = c;
  }

  return HUGE_VAL;
}


/*-----------------------------------------------------------------------------
Name   : newton_raphson()
Purpose: finds a root using Netwon-Raphson method
Returns: the root, or HUGE_VAL if no root was found
Args   : p0       - initial guess
         tol      - tolerance
         max_iter - maximum number of iterations
         f        - the function
         fp       - the derivative of the function
-----------------------------------------------------------------------------*/
double newton_raphson(double p0, double tol, int max_iter,
		      double (*f)(), double(*fp)()) {
  int i = 0;
  double p, divisor;

  printf("Iteration\tApproximation\n"
	 "---------\t-------------\n");

  while (i++ < max_iter) {
    divisor = fp(p0);

    if (divisor == 0.0)
      return HUGE_VAL;

    p = p0 - f(p0) / divisor;

    printf("%9d\t%.12f\n", i, p);

    if (fabs(p - p0) < tol)
      return p;

    p0 = p;
  }

  return HUGE_VAL;
}

/*-----------------------------------------------------------------------------
Name   : secant()
Purpose: finds a root using secant method
Returns: the root, or HUGE_VAL if no root was found
Args   : p0, p1   - initial guesses
         tol      - tolerance
         max_iter - maximum number of iterations
         f        - the function
-----------------------------------------------------------------------------*/
double secant(double p0, double p1, double tol, int max_iter, double (*f)()) {
  int i = 0;
  double p, q0 = f(p0), q1 = f(p1);

  printf("Iteration\tApproximation\n"
	 "---------\t-------------\n");

  while (i++ < max_iter) {

    if (q1 - q0 == 0.0)
      return HUGE_VAL;

    p = p1 - q1 * (p1 - p0) / (q1 - q0);

    printf("%9d\t%.12f\n", i, p);

    if (fabs(p - p1) < tol)
      return p;

    p0 = p1;
    q0 = q1;
    p1 = p;
    q1 = f(p);
  }

  return HUGE_VAL;
}