/* 
 *  Example of using function parameters. 
 * 
 *  Finds roots of the equations 
 *       g(x) = 0   
 *  on a specified interval [x_left, x_right] using the bisection method. 
 *  (cf. Figure 7.10, page 352, HK 5ed.) 
 */ 
 
#include <stdio.h> 
#include <math.h> 

#define FALSE 0
#define TRUE 1

double bisect(double x_left, double x_right, double epsilon,
              double f(double farg), int *errp);
double gcubic(double x);

int 
main(void) 
{ 
    double x_left, x_right, epsilon, root; 
    int    error; 
 
    printf("\nEnter interval endpoints> "); 
    scanf("%lf%lf", &x_left, &x_right); 
    printf("\nEnter tolerance> "); 
    scanf("%lf", &epsilon); 
 
    root = bisect(x_left, x_right, epsilon, gcubic, &error); 

    if (!error) 
        printf("\n gcubic(%.7f) = %e\n", root, gcubic(root)); 

    return (0); 
} 

/* 
 *  Implements the bisection method for finding a root of a function f.
 * 
 *    Bisection method: 
 *      if signs of fp(x_left) and fp(x_right) are different
 *          finds a root and sets output parameter error flag to FALSE 
 *      else 
 *          sets output parameter error flag to TRUE. 
 */ 
double 
bisect(double x_left, double x_right, double epsilon, double f(double farg), 
       int   *errp)       
{ 
    double x_mid,  f_left, f_mid, f_right;  
    int    root_found = FALSE;

    f_left = f(x_left); 
    f_right = f(x_right); 

    if (f_left * f_right > 0) {  
        *errp = TRUE; 
        printf("\nMay be no root in [%.7f, %.7f]", x_left, x_right); 
    } else {   
        *errp = FALSE; 
        while (fabs(x_right - x_left) > epsilon  &&  !root_found) { 

            x_mid = (x_left + x_right) / 2.0;
            f_mid = f(x_mid);
 
            if (f_mid == 0.0)  {          
                root_found = TRUE; 
            } else if (f_left * f_mid < 0.0) {
                x_right = x_mid;                        
            } else {                         
                x_left = x_mid;              
                f_left = f_mid; 
            } 
 
            if (root_found) 
                printf("\nRoot found at x = %.7f, midpoint of [%.7f, %.7f]",
                          x_mid, x_left, x_right); 
            else 
                printf("\nNew interval is [%.7f, %.7f]", x_left, x_right); 
        } 
    } 
 
    return ((x_left + x_right) / 2.0); 
} 
 
/*  
 *  Functions for which roots are sought                             
 */ 
double 
gcubic(double x) 
{ 
   return (5 * pow(x, 3.0) - 2 * pow(x, 2.0) + 3); 
}