//  secant.cpp

   #include <iostream.h>
   #include <iomanip.h>
   #include <math.h>
   #include <string.h> 
   #include <time.h>   

   #define MIN_TOLERANCE 1.0e-10

   int    secant (double, double, int, int &, double &, double &, char
[]); 
   void   get_conditions (double &, double &, int &);
   double user_function (double x);

   int main (void)
   {                               
     clock_t start, finish;
     double  root, value, time_taken, estimate, tolerance;
     int     iterations, max_iterations;
     char    error, err_string[80];
                        
     get_conditions(estimate, tolerance, max_iterations);
                        
     start = clock();

     // Do the calculation 50000 times so that the time taken is a few
     // seconds, rather than a few milliseconds.
     for (unsigned long i = 0 ; i < 50000; i++) 
     {
       error = secant(estimate, tolerance, max_iterations, iterations, 
                      root, value, err_string);
     }            
     
     finish = clock();
     time_taken = (finish - start) / double(CLOCKS_PER_SEC); 

   // Print out the number of iterations, or an error message } 
     
     if (error)
     {
       cout << "ERROR : " << err_string << endl;
     }
     else
     {
       cout << "Successful calculation, " << iterations
            << " iterations." << endl
            << "The root is " << root << endl
            << "The value of the function at this point is "
            << setiosflags(ios::fixed) << setprecision(2) 
            << value << endl
            << "Time for 1000 secant calculations was " 
            << time_taken << " seconds." << endl;         
     }
     
     return 0;
   }


//====================================================================== 

   void get_conditions (double & estimate, double & tolerance, 
          int & max_iterations)
   {     
     do
     {
       cout << "Estimated value of the root (not zero) : ";
       cin  >> estimate; 
     }
     while (estimate == 0.0);
     
     cout << "Calculation will terminate when f(x) is less than the"
          << endl << "tolerance level input by the user." << endl;
     
     do
     {
       cout << "Tolerance                    : ";
       cin  >> tolerance;
       if (tolerance <= MIN_TOLERANCE)
       {
         cout << "ERROR : tolerance must be greater than "
              << MIN_TOLERANCE << endl;
       }
     }
     while (tolerance <= MIN_TOLERANCE);

     do
     {
       cout << "Maximum number of iterations (1 or more) : ";
       cin  >> max_iterations;
     }
     while (max_iterations < 1);
   }
          


//====================================================================== 
   // Function to find a root of an equation by the secant method. 
   // This is basically the same as the Newton method, but the 
   // approximation is made that dy/dx = (x2-x1)/(y2-y1), so the der
   // is not directly calculated.
   //   estimate       : Initial estimate for the root.  This is x1 in the 
   //                    approximation above - x2 is the estimate plus 1%
   //   tolerance      : tolerance level of the answer
   //   max_iterations : maximum number of iterations to be performed
   //   root           : final value of the root
   //   value          : value of user_funct at the root
   //   iterations     : number of iterations performed 
   //   err_string     : string containing an error message, if approp.
   // Function returns 1 (true) if an error has been encountered. 

//====================================================================== 

   int  secant (double estimate, double tolerance, int max_iterations, 
          int & iterations, double & root, double & value, char
err_string[])
   { 
     double y1, y2, y3, a1, a2, a3;
     
     a1 = estimate;         y1 = user_function(a1);
     a2 = estimate * 1.01;  y2 = user_function(a2);

   // Begin the iterations.

     for (iterations = 1; iterations <= max_iterations; iterations++)
     {
       if (fabs(y1 - y2) < 1.0e-10)
       {
         strcpy(err_string, "y1 = y2");
         return 1;
       }
       
       a3 = a2 - y2*(a2 - a1)/(y2 - y1);
       y3 = user_function(a3);
       a1 = a2;  y1 = y2;
       a2 = a3;  y2 = y3;
       
       if (fabs(y3) <= tolerance)
       {
         root  = a3;
         value = y3;
         return 0;
       }
     }
     
     strcpy(err_string, "Maximum number of iterations was exceeded");
     return 1;
   }


//====================================================================== 
   // This function describes a mathematical function whose root is
   // to be found.

//====================================================================== 

   double user_function (double x)
   {
     return cos(x) - x;
   }