/********
 * find_zero.c
 * 
 *  finds the zero of a function
 *  by bisecting the range where
 *  the zero is to be found on 
 *  each iteration.
 *
 *  compile with "cc -lm find_zero_FIRST.c -o find_zero_FIRST"
 * 
 *  ===========================================================================
 *  == This program has a (at least one) serious bug in it - can you find it? =
 *  ===========================================================================
 *  
 * 
 *   - Jim Mahoney, Sep 2001
 ***********/       
#include <stdio.h>
#include <math.h>
#include <float.h>


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


main () {
  double B, accuracy;
  double f1, f2, fmid;
  double x1, x2, xmid;
  int iteration = 0;
  int iterationMax = 100; /* exit if we loop more than this */

  printf("\n -- find_zero.c --  \n\n");

  /* -- Read in the inputs -- */

  printf("What is the parameter B ? \n");	scanf("%lf", &B);
  printf("What is xMin ? \n");			scanf("%lf", &x1);
  printf("What is xMax ? \n");			scanf("%lf", &x2);
  printf("What accuracy do you want? \n");	scanf("%lf", &accuracy);

  printf("  (B, xMin, xMax, accuracy) = (%g,%g,%g,%g)\n", B,x1,x2,accuracy);

  f1 = f(x1,B);
  f2 = f(x2,B);

  /* -- Consistency checks -- */

  if ( f1*f2 > 0.0 ) {
    printf(" ** ERROR ** \n");
    printf(" f(xMin,B) and f(xMax,B) have the same sign. \n");
    exit();
  }

  accuracy = abs(accuracy);	/* make sure accuracy is positive */

  /* -- loop until the endpoints are close enough -- */

  while ( abs(x1-x2) > accuracy ) {
    iteration++;

    if (iteration>100){
      printf(" ** ERROR ** \n  too many iterations \n");
      exit();
    }

    /* The heart of the calculation. */
    /* Do you see how it works?      */
    xmid = (x1+x2)/2.0;
    fmid = f(xmid,B);
    if (fmid*f1 > 0) {
      x1 = xmid;
      f1 = fmid;
    }
    else {
      x2 = xmid;
      f2 = fmid;
    }
  }

  printf(" Done. \n");
  printf("  f( %.20g, %g ) = %g \n", x1, B, f1);
  printf("  f( %.20g, %g ) = %g \n", x2, B, f2);
  printf("\n");
  printf(" best guess midpoint = %.20g \n", (x1+x2)/2.0 );

  printf("\n\n");
  printf(" -- from float.h -- \n");
  printf(" DBL_EPSILON = smallest such that 1.0 +x != 1.0 : %g \n", DBL_EPSILON);
  printf(" DBL_DIG = decimal digits of precision: %i \n", DBL_DIG);

}