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

/* signum: returns the sign of x; returns 0 if x.eq.0. */
double signum(double x)
{
  if (x<0) {
    return -1;
  }
  else if (x==0) {
    return 0;
  }
  else if (x>0) {
    return 1;
  }
  else {
    /* must be that x is Not-a-Number */
    return x;
  }
}

/* func0: test function for rootfinding routine. */
/* func0 is smooth and has a simple root at x=sqrt(2.0). */
double func0(double x)
{
  return (x-sqrt(2.0))*sqrt(1+x*x);
}

/* func1: test function for rootfinding routine. */
/* func1 is smooth and has a triple root at x=sqrt(2.0). */
double func1(double x)
{
  return pow((x-sqrt(2.0)),3);
}

/* func2: test function for rootfinding routine. */
/* func2 is smooth and has a simple root at x=sqrt(2.0); the */
/* asymptotic behavior may pose problems for a Newton-like method. */
double func2(double x)
{
  return atan((x-sqrt(2.0)));
}

/* func3: test function for rootfinding routine. */
/* func3 is smooth and has a simple root at x=sqrt(2.0); the */
/* asymptotic behavior may pose problems for a Newton-like method. */
double func3(double x)
{
  return atan(256*(x-sqrt(2.0)));
}

/* func4: test function for rootfinding routine. */
/* func4 is discontinuous with sign change at x=sqrt(2.0). */
double func4(double x)
{
  if (x<0||x*x<2) {
    return -1;
  }
  else {
    return 1;
  }
}

/* func5: test function for rootfinding routine. */
/* func5 has a vertical asymptote and sign change at x=sqrt(2.0). */
double func5(double x)
{
  double t;
  t = (x-sqrt(2.0))*sqrt(1+x*x);
  if (t!=0) {
    return 1/t;
  }
  else {
    return 0;
  }
}

/* func6: test function for rootfinding routine. */
/* func6 is discontinous with sign change at 0. */
double func6(double x)
{
  if (x<0) {
    return -1;
  }
  else {
    return 1;
  }
}

/* findroot: locate a root or a sign change of a given function in
/* a given interval.
/* Inputs:
/*   f: double function of one double argument.
/*   a0, b0: double, bounding bracket for the root.
/*   Must have signum(f(a0))*signum(f(b0))<=0.
/* Returns: 
/*   xend: root of f in the interval defined by (a0,b0).
/*   nmbf: number of calls to f made for this call to findroot.
*/
void findroot (double (*f)(double), double a0, double b0,
	       double *xend, int *nmbf)
{
  ...
}

main()
{
  double a0, b0, xend;
  int nmbf;
  printf("  %8s  %16s  %16s  %6s\n", "function", "xend", "error", "nmbf");
  a0=1; b0=2;
  findroot(&func0, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func0", xend, xend-sqrt(2.0), nmbf);
  a0=1; b0=2;
  findroot(&func1, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func1", xend, xend-sqrt(2.0), nmbf);
  a0=1; b0=2;
  findroot(&func2, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func2", xend, xend-sqrt(2.0), nmbf);
  a0=1; b0=2;
  findroot(&func3, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func3", xend, xend-sqrt(2.0), nmbf);
  a0=1; b0=2;
  findroot(&func4, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func4", xend, xend-sqrt(2.0), nmbf);
  a0=1; b0=2;
  findroot(&func5, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func5", xend, xend-sqrt(2.0), nmbf);
  a0=-1; b0=2;
  findroot(&func6, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func6", xend, xend-0, nmbf);
}