/* Find a root of f() by the bisection method. */

#define   EPS   1e-12               /* epsilon */

double root(double f(double), double a, double b)
{
   double   m = (a + b) / 2.0;     /* midpoint */

   if (f(m) == 0.0 || b - a < EPS)
      return m;
   else if (f(a) * f(m) < 0.0)
      return root(f, a, m);
   else
      return root(f, m, b);
}