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

/* signum: returns the sign of x; returns 0 if x==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;
  }
}

/* bisect: returns a bisection point for an interval.
   If a and b are both non-zero, of the same sign, and not too far
   apart on a logarithmic scale, returns c=(a+b)/2.  Other cases ...
   see the code.  The intent is that bisect(a,b) divides the interval
   so that there are about equal number of floating point values in
   either half.  Some preference towards zero. */
double bisect(double a, double b)
{
  double eps, tiny, t0, t1;
  eps = pow(2.0,-52); /* IEEE machine precision */
  tiny = pow(2.0,-1022); /* IEEE smallest positive normal number */
  if (signum(a)*signum(b)<0) {
    return 0;
  }
  else if (signum(a)*signum(b)>0) {
    if (fabs(a)<fabs(b)) {
      t0 = fabs(a); t1 = fabs(b);
    }
    else {
      t0 = fabs(b); t1 = fabs(a);
    }
    if (t1<4*t0||t1<4*tiny) {
      return (a+b)/2;
    }
    else {
      return signum(a)*sqrt(t0)*sqrt(t1);
    }
  }
  else if (a==0||b==0) {
    t0 = a+b;
    if (fabs(t0)>=1/sqrt(eps)) {
      return signum(t0);
    }
    else if (fabs(t0)>=sqrt(eps)) {
      return t0*sqrt(eps);
    }
    else if (fabs(t0)>=sqrt(tiny)) {
      return t0*fabs(t0);
    }
    else if (fabs(t0)>=4*tiny) {
      return signum(t0)*tiny;
    }
    else if (fabs(t0)>eps*tiny) {
      return signum(t0)*eps*tiny;
    }
    else {
      return 0;
    }
  }
  else {
    /* must be we have a NaN */
    return (a+b)/2;
  }
}

/* func0: test function for rootfinding routine.
   func0 is smooth and has a simple root at x=sqrt(2.0d0). */
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.0d0). */
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.0d0); 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.0d0); 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.0d0). */
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.0d0). */
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_bisect: locate a root or a sign change of a given function in
   a given interval. */
void findroot_bisect (double (*f)(double), double a0, double b0,
	       double *xend, int *nmbf)
{
  double a, b, c, fa, fb, fc;
  a = a0; fa = f(a);
  b = b0; fb = f(b);
  *nmbf = 2;
  if (fa==0) {
    c = a; fc = fa;
  }
  else if (fb==0) {
    c = b; fc = fb;
  }
  else {
    assert(signum(fa)!=signum(fb));
    c = bisect(a,b);
    fc = f(c);
    ++*nmbf;
  }
  while (c!=a&&c!=b) {
    if (fc==0) {
      a = c; fa = fc;
      b = c; fb = fc;
    }
    else if (signum(fc)==signum(fa)) {
      a = c; fa = fc;
      c = bisect(a,b);
      fc = f(c);
      ++*nmbf;
    }
    else if (signum(fc)==signum(fb)) {
      b = c; fb = fc;
      c = bisect(a,b);
      fc = f(c);
      ++*nmbf;
    }
    else {
      assert(1);
    }
  }
  *xend = c;
  return;
}

main()
{
  double a0, b0, xend;
  int nmbf;
  printf("  Bisection method\n");
  printf("  %8s  %16s  %16s  %6s\n", "function", "xend", "error", "nmbf");
  a0=1; b0=2;
  findroot_bisect(&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_bisect(&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_bisect(&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_bisect(&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_bisect(&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_bisect(&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_bisect(&func6, a0, b0, &xend, &nmbf);
  printf("  %8s  %16.7e  %16.7e  %6d\n", "func6", xend, xend-0, nmbf);
}