#include
#include


/* Program to demonstrate numerical root finding methods.  Starting
  from some initial guesses at a root (i.e. value of x such that
  F(x) = 0), these algorithms iteratively improve the estimate of
  the root until some convergence criterion is fulfilled.

  Three techniques are presented.  The Bisection method is easy to
  understand and safe, but converges slowly.  The Newton-Raphson
  method converges quickly, but requires evaluation of the derivative
  of F.  The Secant method converges nearly as quickly as Newton-
  Raphson, but does not require any derivative evaluations.
*/

#define MinTol 1.0E-10
#define TrueRoot 0.772883
#define MaxStep 20



double f(double x){
/* The function for which a root is to be found */
  return(x*x*x - exp(-x));
}


double f_prime(double x){
/* The derivative of f */
  return(3*x*x + exp(-x));
}


void swap(double *a, double *b){
/* Exchange two real variables */

  double temp;
  temp = *a;
  *a = *b;
  *b = temp;
}



double bisect(double x0, double x1, double tolerance){

/* Finds one real root of the function F(x) on the interval
  [X_Lo, X_Hi], which must contain a sign change in F (i.e.
  this will not work for roots of even multiplicity).

  Proceeds by halving the interval [X_Lo, X_Hi].  Depending
  on the sign of the function at the midpoint of the interval,
  either the low or high limit is discarded (in order to preserve
  the sign change in the interval).

  The algorithm terminates when either the size of the interval
  shrinks below Tolerance, or the magnitude of F at the midpoint
  of the interval is less than Tolerance. */

  double root;
  int i;

  for(i=1; i<=MaxStep; i++){
      root=(x0 + x1) /2;
      if(f(x0)*f(root) > 0.0)
        x0=root;
      else
        x1=root;
      if( (x1-root) < tolerance || fabs(f(root)) < tolerance) break;
      }
  return(root);
}


double newton_raphson(double x, double tolerance){
/* Finds one real root of the function F(x) using the Newton-
  Raphson method and an initial guess at the root, x0.

  While this method converges rapidly to the root, it does
  require the evaluation of the derivative of F, F'(x).

  The algorithm terminates when either the difference between
  new and old estimates of the root falls below Tolerance, or
  the magnitude of the function at the latest estimate is less
  than Tolerance.

  As a fail-safe, the algorithm is abandoned after MaxStep
  (a global constant) iterations. */

  double root;
  int i;

  for(i=1; i<=MaxStep; i++){
      root = x - (f(x) / f_prime(x));
      if( (root-x) < tolerance || fabs(f(root)) < tolerance) break;
      x=root;

      }
  return(root);
}


double secant(double x0, double x1, double tolerance){
/* Finds one real root of the function F(x) using the Secant
  method and two initial guesses at the root, x0 and x1.

  This method has similar convergence properties to Newton-
  Raphson, but does not require the evaluation of F'.

  The algorithm terminates when either the difference between
  new and old estimates of the root falls below Tolerance, or
  the magnitude of the function at the latest estimate is less
  than Tolerance.

  As a fail-safe, the algorithm is abandoned after MaxStep
  (a global constant) iterations.

  If you understand the Newton-Raphson method, can you see how
  the f(x0) * (x2-x1)/(f(x2) - f(x1)) term in this method replaces
  the f(x)/f_prime(x) ?  This secant method is actually approximating
  the slope, to avoid evaluating the derivative (which may at times
  be impossible to do.) */


  double root;
  int i;

  for(i=1; i<=MaxStep; i++){
      root = x1 - f(x0)*(x1-x0)/(f(x1)-f(x0));
      if( fabs(x1-x0) < tolerance || fabs(f(root)) < tolerance) break;
      x0 = x1;

      x1= root;

      }
  return(root);
}

main(){

  double x0, x1, tolerance, root;

  printf("Root-finding Program\n\n\n");
  printf(" This program finds the root of the function x^3 - e^(-x).  It\n");
  printf("will prompt you for two guesses, as well as the tolerance, and\n");
  printf("then it will use the bisection method, the secant method, and \n");
  printf("the Newton-Raphson method to find the root.  The actual root \n");
  printf("is approximately  x = 0.772883  .\n\n\n");
  printf("Enter first point value  x=");
  scanf("%lf",&x0);
  printf("Enter last point value  x=");
  scanf("%lf",&x1);
  printf("Enter Tolerance :");
  scanf("%lf",&tolerance);

  if(x0 > x1) swap(&x0, &x1);

  if( fabs(f(x0)) <= tolerance || fabs(f(x1)) <= tolerance)

      printf("One of the endpoints is a root.\n");

  else if( f(x0)*f(x1) > 0 )
      printf("There is no single root in that interval.\n");

  else if(tolerance < MinTol)

      printf("That is too small a tolerance.\n");

  else if( (x1 - x0) < tolerance )
      printf("The specified interval is smaller than the tolerance.\n");

  else{

  /* Approximate the root with bisection method */
      root=bisect(x0, x1, tolerance);
  /* If the root returned does not yield zero within tolerance,
     it means the function could not find the root with the allowed
     number of iterations. */
      if(f(root) > tolerance){
        printf("The bisection method could not find a root within\n");
        printf("tolerance with %d iterations.\n",MaxStep); }
      else{
        printf("\nBisection method yields;\n");
        printf("Root=%8.5lf True Root=%8.5lf Error=%8.5lf",
               root, TrueRoot, (root-TrueRoot) );  }

  /* Approximate the root with Newton-Raphson method */
      root=newton_raphson(x0, tolerance);
      if(f(root) > tolerance){
        printf("\nThe Newton-Raphson method could not find a root within\n");
        printf("tolerance with %d iterations.\n",MaxStep); }
      else{
        printf("\nNewton-Raphson method yields;\n");
        printf("Root=%8.5lf True Root=%8.5lf Error=%8.5lf",
               root, TrueRoot, (root-TrueRoot) );}

  /* Approximate the root with the secant method */
      root=secant(x0, x1, tolerance);
      if(f(root) > tolerance){
        printf("\nThe secant method could not find a root within\n");
        printf("tolerance with %d iterations.\n",MaxStep); }
      else{
        printf("\nSecant method yields;\n");
        printf("Root=%8.5lf True Root=%8.5lf Error=%8.5lf\n",
               root, TrueRoot, (root-TrueRoot) );  }

  } /* END the long if-else block */

}  /* END main() */