/********
 * zeros
 * 
 *  finds the zero of a function
 *  by bisecting the range where
 *  the zero is to be found on 
 *  each iteration.
 *
 *  compile with "cc -lm zeros.c -o zeros"
 * 
 *   - Jim Mahoney, Sep 2001
 ***********/       
#include <stdio.h>
#include <math.h>

/****************************************************
/* Here's the function that we'll find the zero of. */
double f(double x, double B) {
  return( sin(x) - B*x*x );
}

/* and its derivative  */
double dfdx(double x, double B) {
  return( cos(x) - 2.0*B*x );
}

/* and a way to print it out for the user to see. */
void print_f(){
  printf(" sin(x) - B*x*x " );
}

/*************************************************
/* Find the value of x1< x < x2 where f(x,B) = 0 */
/* using a bisection algorithm, continuing until the */
/* error in x stops decreasing or we hit the iteration limit. */
double bisection( double B, double x1, double x2) {
  double f1,f2,dx,dxLast,xmid,fmid;
  int iterationMax = 100;
  int i=0;
  f1 = f(x1,B);
  f2 = f(x2,B);
  if ( f1*f2 > 0.0 ) {    /* -- Consistency checks -- */
    printf(" ** ERROR in bisection ** \n");
    printf(" f(xMin,B) and f(xMax,B) have the same sign. \n");
    exit();
  }
  dx = fabs(x1-x2);  
  dxLast = 2.0*dx;
  while ( (dx < dxLast) && (i<iterationMax) ) {
    i++;
    printf(" i=%i: (x1,x2)=(%g,%g); (f1,f2)=(%g,%g) \n", i,x1,x2,f1,f2);
    dxLast = dx;
    xmid = (x1+x2)/2.0;
    fmid = f(xmid,B);
    if (fmid*f1 > 0) {
      x1 = xmid;
      f1 = fmid;
    }
    else {
      x2 = xmid;
      f2 = fmid;
    }
    dx = fabs(x1-x2);
  }
  return( (x1+x2)/2.0 );
}

/*************************************************
/* Find the zero of f(x,B) using its first derivative dfdx
/* and a starting point x0 , continuing until the */
/* error in x stops decreasing or we hit the iteration limit. */
double newtons( double B, double x0) {
  double fx0, deriv, dx, dxLast, xNew;
  int iterationMax = 100;
  int i=0;
  dx = 1e20;      /* first step should be smaller than this */
  dxLast =  2e20; /* which should be smaller than this */
  while ( (dx < dxLast) && (i<iterationMax) ) {
    i++;
    dxLast = dx;
    deriv = dfdx(x0,B);
    fx0 = f(x0,B);
    printf(" i=%i: (x,f,deriv)=(%.16g,%g,%g) \n", i, x0, fx0, deriv );
    if ( deriv == 0.0 ) {
      printf( " error : deriv=0 in newtons \n ");
      exit();
    }
    xNew = x0 - fx0 / deriv;
    dx = fabs(x0-xNew);
    x0 = xNew;
  }
  return( x0 );
}




/* --------------------------------------------------------------- */
main () {
  double B, x, fx;
  double xmin, xmax, x0;
  double Bmin, Bmax, dB;  /* */

  printf("\n -- zeros --  \n\n");

  /* -- set the parameter -- */
  xmin = 1e-6;  /* just bigger than 0 */
  xmax = 3.14;  /* just less than pi */
  x0 = 1.0;     /* something near the zero we want */
  B = 1.0;

  x = bisection(B,xmin,xmax);
  printf( " bisection : f(x=%.16g,B=%g) = %g \n\n", x, B, f(x,B) );

  x = newtons(B,x0);
  printf( " newtons : f(x=%.16g,B=%g) = %g \n\n", x, B, f(x,B) );



}