////////////////////////////////////////////////////////////////////////////////
// File: secant_method.c                                                      //
// Routines:                                                                  //
//    double Secant_Method                                                    //
////////////////////////////////////////////////////////////////////////////////

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

////////////////////////////////////////////////////////////////////////////////
//  double Secant_Method( double (*f)(double), double a, double b,            //
//                      double tolerance, int max_iteration_count, int *err)  //
//                                                                            //
//  Description:                                                              //
//     Estimate the root (zero) of f(x) using the secant method where 'a'     //
//     and 'b' are two initial estimates.  The process terminates successfully//
//     when the change in the estimate of the root is less than 'tolerance'   //
//     or unsuccessfully when the number of iterations reaches 'max_iteration_//
//     count' or if an attempt is made to divide by zero.                     //
//                                                                            //
//     The secant method uses two points on the curve f(x) to estimate the    //
//     root as the intersection of the x-axis with the line joining the       //
//     two points (the secant).  The process then replaces the oldest of the  //
//     two estimates with the newest estimate.  The process continues until   //
//     the distance between the two estimates is less than the tolerance.     //
//     The process can fail if there is a local maximum or minimum near the   //
//     the root and is sensitive to having good initial estimates.  If there  //
//     is a local maximum or minimum in the neigborhood of the root, the      //
//     return value should be verified before use.  The secant method is      //
//     similar to the Newton-Raphson method in which the derivative is        //
//     replaced by a numerical estimate of the derviative.                    //
//                                                                            //
//  Arguments:                                                                //
//     double *f         Pointer to function of a single variable of type     //
//                       double.                                              //
//     double  a         Initial estimate.                                    //
//     double  b         Initial estimate.                                    //
//     double  tolerance Desired accuracy of the zero, tolerance > 0.         //
//     int     max_iteration_count Maximum allowable number of iterations.    //
//     int     *err      0 if successful, -1 if not, i.e. the iteration       //
//                       count meets or exceeds the max_iteration_count,      //
//                       -2 if an attempt is made to divide by zero.          //
//                                                                            //
//  Return Values:                                                            //
//     A zero somewhere in a neighborhood of a and b.                         //
//                                                                            //
//  Example:                                                                  //
//     {                                                                      //
//        double f(double), zero, a, b, tolerance = 1.e-6;                    //
//        int max_iteration = 10;                                             //
//        int err;                                                            //
//                                                                            //
//        (determine nearby values a and b of a zero)                         //
//        zero = Secant_Method( f, a, b, tolerance, max_iterations, &err);    //
//        ...                                                                 //
//     }                                                                      //
//     double f(double x) { define f }                                        //
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
double Secant_Method( double (*f)(double), double a, double b, 
                        double tolerance, int max_iteration_count, int *err) {

   double fa = (*f)(a), fb = (*f)(b), delta;
   int i;

   *err = 0;
   for (i = 0; i < max_iteration_count; i++) {

         // Insure that the f(b) has smaller magnitude than f(a). //

      if ( fabs(fa) < fabs(fb) ) {
         delta = a; a = b; b = delta; delta = fa; fa = fb; fb = delta;
      }
    
         // If an attempt to divide by zero, return with an error code -2. //

      if ( fb == fa ) { *err = -2; return 0.0; }

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

      delta = fb * (a - b) / (fb - fa);

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

      if ( delta == 0.0 ) return b;

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

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

      a = b + delta;
      fa = (*f)(a);
   }
   
   *err = -1;
   return a;   
}