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

#include <stdio.h>

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

double p(double);
double root(double (*)(double), double, double);

int main(void)
{
  double   x;

  x = root(p, 0.0, 3.0);
  printf("%s%.16f\n%s%.16f\n",
         "Approximate root: ", x,
         "  Function value: ", p(x));
  return 0;
}

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

/* a polynomial p */
double p(double x)
{
  return (x * x * x * x * x  -  7.0 * x  -  3.0);
}