#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;
  }
}

/* min: returns the minimum of x and y. */
double min(double x, double y)
{
  if (x<=y) {
    return x;
  }
  else if (y<=x) {
    return y;
  }
  else {
    /* must have a NaN */
    return (x+y)/2;
  }
}

/* max: returns the maximum of x and y. */
double max(double x, double y)
{
  if (x<=y) {
    return y;
  }
  else if (y<=x) {
    return x;
  }
  else {
    /* must have a NaN */
    return (x+y)/2;
  }
}

/* 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;
  }
}

/* fraclin: fractional linear interpolation. */
double fraclin(double x0, double y0,
               double x1, double y1,
               double x2, double y2)
{
  double dx0, dx1, xnew0, xnew1, t0, t1;
  dx0 = x0-x2; dx1 = x1-x2;
  xnew0 = x2+dx0*dx1*y2*(y1-y0)/(dx1*y1*(y2-y0)-dx0*y0*(y2-y1));
  xnew1 = x2+dx0*dx1*(y1-y0)/(dx1*(y2-y0)-dx0*(y2-y1));
  t0 = min(min(x0,x1),x2);
  t1 = max(max(x0,x1),x2);
  if (t0<=xnew0&&xnew0<=t1) {
    return xnew0;
  }
  else if (t0<=xnew1&&xnew1<=t1) {
    return xnew1;
  }
  else {
    /* this can only occur because of roundoff error.  Use the median. */
    return min(min(max(x0,x1),max(x1,x2)),max(x2,x0));
  }
}

/* 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;
}

/* findroot: locate a root or a sign change of a given function in
   a given interval. */
void findroot (double (*f)(double), double a0, double b0,
	       double *xend, int *nmbf)
{
  double a, b, c, fa, fb, fc, t0, xnew;
  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)) {
      t0 = fraclin(a,fa,b,fb,c,fc);
      if (fabs(b-c)<2*fabs(b-a)/3&&
	  min(c,b)<t0&&t0<max(c,b)) {
	xnew = t0;
      }
      else {
	xnew = bisect(c,b);
      }
      a = c; fa = fc;
      c = xnew;
      fc = f(c);
      ++*nmbf;
    }
    else if (signum(fc)==signum(fb)) {
      t0 = fraclin(a,fa,b,fb,c,fc);
      if (fabs(c-a)<2*fabs(b-a)/3&&
	  min(a,c)<t0&&t0<max(a,c)) {
	xnew = t0;
      }
      else {
	xnew = bisect(a,c);
      }
      b = c; fb = fc;
      c = xnew;
      fc = f(c);
      ++*nmbf;
    }
    else {
      assert(1);
    }
  }
  *xend = c;
  return;
}

/* 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;
  }
}

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);
  printf("  Fractional linear method\n");
  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);
}