////////////////////////////////////////////////////////////////////////////////
// File: newton_raphson.c                                                     //
// Routines:                                                                  //
//    double Newton_Raphson                                                   //
////////////////////////////////////////////////////////////////////////////////

#include <math.h>                                        // required for fabs()

////////////////////////////////////////////////////////////////////////////////
//  double Newton_Raphson( double (*f)(double), double (*df)(double),         //
//            double a, double tolerance, int max_iteration_count, int *err)  //
//                                                                            //
//  Description:                                                              //
//     Estimate the root (zero) of f(x) using the Newton-Raphson method       //
//     with 'a' as an initial estimate.  The process terminates successfully  //
//     when either the function vanishes or when the change in the estimate   //
//     of the root is less than 'tolerance' or unsuccessfully when the        //
//     number of iterations reaches 'max_iteration_count'.                    //
//                                                                            //
//     The Newton-Raphson method uses the intersection of the x-axis and the  //
//     tangent to the curve at an estimate of the root to obtain a new        //
//     estimate of the location of the root.  The process is continued until  //
//     the change of locations of the estimates is within a preassigned       //
//     tolerance.  The Newton-Raphson method requires the user not only       //
//     furnish the function but also its derivative.  The method can fail if  //
//     the root is near a local minimum or a local maximum.  If the order     //
//     of the zero of f(x) is greater than one, one can find a root of        //
//     f(x) / f'(x).  The Newton-Raphson algorithm is quadratically           //
//     convergent near a simple root.                                         //
//                                                                            //
//  Arguments:                                                                //
//     double *f         Pointer to function of a single variable of type     //
//                       double.                                              //
//     double *df        Pointer to the derivative of the function *f.        //
//     double  a         Initial estimate.                                    //
//     double  tolerance Desired accuracy of the zero.                        //
//     int     max_iteration_count Terminate iteration if the number of       //
//                                 iterations meets or exceed this number.    //
//     int     *err      0 if successful, -1 if not, i.e. the number of       //
//                       iterations meets or exceeds the max_iteration_count, //
//                       or -2 if the derivative vanishes at an evaluation    //
//                       point.                                               //
//                                                                            //
//  Return Values:                                                            //
//     A root somewhere in a neighborhood of a.                               //
//                                                                            //
//  Example:                                                                  //
//     {                                                                      //
//        double f(double), df(double),zero, a, tolerance = 1.e-6;            //
//        int max_iteration = 10;                                             //
//        int err;                                                            //
//                                                                            //
//        (determine the initial estimate a of the root.)                     //
//        zero = Newton_Raphson( f, df, a, tolerance, max_iteration, &err);   //
//        ...                                                                 //
//     }                                                                      //
//     double f(double x) { define f }                                        //
//     double df(double x) { define df }                                      //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
double Newton_Raphson( double (*f)(double), double (*df)(double), double a,
                       double tolerance, int max_iteration_count, int *err) {

   double delta, fa, dfa;
   int i;
   
   *err = 0;
   for (i = 0; i < max_iteration_count; i++) {
      if ( ( fa = (*f)(a) ) == 0.0 ) return a;

         // If an attempt to divide by zero, return with an error code -2. //

      if ( ( dfa = (*df)(a) ) == 0.0 ) { *err = -2; return 0.0; }

         // Estimate the location of the root relative to a. //

      delta = - fa / dfa;

         // If the location of the root relative to a is less than //
         // the tolerance, return the estimate of the root.        //

      if (fabs(delta) < tolerance ) return a + delta;

         // Otherwise, update the estimate of the root. //

      a += delta;
   }

   *err = -1;
   return a;
}