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

#define EPSILON 1E-9


double
eval_poly(double x, double coeffs[], int order)
{
  double res;
  int i;

  res = coeffs[0];

  for ( i = 1; i <= order; i++ ) {
    res = res * x + coeffs[i];
  }

  return res;
}


double
find_zero(double left, double right, double coeffs[], int order)
{
  double c, lp, rp, cp, eps=EPSILON;

  if(eval_poly(left, coeffs, order) * eval_poly(right, coeffs, order) >= 0) {
    exit(1);
  }

  while(1) {
    c = (left+right)/2;

    lp = eval_poly(left, coeffs, order);
    rp = eval_poly(right, coeffs, order);
    cp = eval_poly(c, coeffs, order);

    if(fabs(cp) < eps) {
      break;
    }

    if(cp * lp < 0) {
      right = c;
    }
    else if(cp * rp < 0) {
      left = c;
    }
    else {
      exit(1);
    }
  }

  return c;
}

int
main(void)
{
  double coeffs[] = { 7, 5, 3, 2, 1, 1 };

  printf("eval_poly: %f\n", eval_poly(3, coeffs, 5));
  printf("%f\n", find_zero(-100, 100, coeffs, 5));

  return 0;
}