(wrap:c-lines "
#include <math.h>
#include <stdio.h>

#define SIGNUM(x) (((x) > 0) ? 1 : ((x) < 0) ? -1 : 0)

double zero(double (*f)(double), double a, double b, double fa, double fb,
            double tol)
{
  if (b - a <= tol)
    return (b + a)/2.0;
  else if (0 <= fa && 0 <= fb)
    return fa < fb ? a : b;
  else if (0 >= fa && 0 >= fb)
    return fa < fb ? b : a;
  else {
    double c = (b + a) / 2.0;
    double fc = f(c);
    if (SIGNUM(fa) == SIGNUM(fc))
      return zero(f, c, b, fc, fb, tol);
    else
      return zero(f, a, c, fa, fc, tol);
  }
}
")
(wrap:c-callback "g" zerofun (:flonum) :flonum)
(wrap:c-lines "
double base_zero(double a, double b, double tol)
{
  return zero(g, a, b, g(a), g(b), tol);
}
")
(wrap:c-function base-zero "base_zero" (:flonum :flonum :flonum) :flonum)
(defvar *zerofun*)
(defun zerofun (x) (funcall *zerofun* x))
(defun zero (f a b &optional (tol .00001))
  (let ((*zerofun* f))
    (base-zero a b tol)))
(wrap:c-callback-variable "gv" *zerofun* (:flonum) :flonum)
(wrap:c-lines "
double base_zero_v(double a, double b, double tol)
{
  return zero(gv, a, b, gv(a), gv(b), tol);
}
")
(wrap:c-function base-zero-v "base_zero_v" (:flonum :flonum :flonum) :flonum)
(defun zero-v (f a b &optional (tol .00001))
  (let ((*zerofun* f))
    (base-zero-v a b tol)))