#include <stdio.h>
#include <math.h>

/*
Compare the performance of bisection, the secant method, and
Newton's method for root finding.
*/

/* The objective function and derivative. The root is x=pi. */
double fun(const double x) { return tan(x/4) - 1; }
double der(const double x) { return 1 / cos(x/4) / cos(x/4) / 4; }

/* The sign function. */
int sgn(const double x) { return (x > 0) ? 1 : -1; }


/* Take one step of the bisection algorithm.  The bracket abscissas
   are x1,x2 and the bracket ordinates are y1,y2. */
void bisect_step(double* x1,
		 double* x2,
		 double* y1,
		 double* y2)
{
  double x = (*x1+*x2)/2;
  double y = fun(x);
  if (sgn(y) == sgn(*y1)) { *x1 = x; *y1 = y; }
  else { *x2 = x; *y2 = y; }
}


/* Make one step of the secant method.  The decant is drawn between x1
 and x2.  On entry, y1 = f(x1) and y2 = f(x2).  On exit, either
 (x1,y1) or (x2,y2) is updated to hold the new endpoint. */
void secant_step(double* x1,
		 double* x2,
		 double* y1,
		 double* y2)
{
  double d = *y1-*y2, x;

  if (fabs(d) > 1e-16)
    {
      x = *x1 + *y1*(*x2-*x1)/d;
      *x2 = *x1;
      *y2 = *y1;
      *x1 = x;
      *y1 = fun(*x1);
    }
}    


int main(int argc, char** argv)
{
  int i;
  double E, E1, x1, x2, y1, y2, nx;
  const double PI = 4*atan(1);

  /*
    Test bisection.
  */

  /* A bracket. */
  x1=0; 
  x2=4; 
  y1=fun(x1); 
  y2=fun(x2);
  
  /* Run 50 steps of bisection. */
  printf("50 steps of bisection:\n");
  printf("Value                    Error                    "
	 "Error/last error\n");
  for (i=0; i<50; ++i)
    {
      bisect_step(&x1, &x2, &y1, &y2);
      E = PI - x1;
      printf("%20.17f\t%20.17f\t%20.16f\n", x1, E, (i == 0) ? 0 : E/E1);
      E1 = E;
    }
  printf("\n\n");

  /*
    Test the secant method.
  */

  /* A bracket. */
  x1 = 0;
  x2 = 4;
  y1 = fun(x1); 
  y2 = fun(x2);

  /* Run 50 steps of the secant method. */
  printf("50 steps of the secant method:\n");
  printf("Value                    Error                    "
	 "Error/last error\n");
  for (i=0; i<50; ++i)
    {
      secant_step(&x1, &x2, &y1, &y2);
      E = PI - x1;
      printf("%20.17f\t%20.17f\t%20.16f\n", x1, E, (i == 0) ? 0 : E/E1);
      E1 = E;
    }
  printf("\n\n");

  /*
    Test Newton's method.
  */
  
  /* Starting point for Newton's method. */
  nx = 0;
  
  /* Run 50 steps of Newton's method. */
  printf("50 steps of Newton's method:\n");
  printf("Value                    Error                    "
	 "Error/last error\n");
  for (i=0; i<50; ++i)
    {
      nx -= fun(nx) / der(nx);
      E = PI - nx;
      printf("%20.17f\t%20.17f\t%20.16f\n", nx, E, (i == 0) ? 0 : E/E1);
      E1 = E;
    }
  printf("\n\n");

  return 0;
}