/* zeroin.f -- translated by f2c (version 20000817).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

doublereal zeroin_(doublereal *ax, doublereal *bx, doublereal (*f) (), doublereal *tol, integer flag)

{
    /* System generated locals */
    doublereal ret_val, d__1;

    /* Builtin functions */
    double d_sign();

    /* Local variables */
    static doublereal a, b, c__, d__, e, p, q, r__, s, fa, fb, fc, xm, eps, 
	    tol1;


/*      a zero of the function  f(x)  is computed in the interval ax,bx . */

/*  input.. */

/*  ax     left endpoint of initial interval */
/*  bx     right endpoint of initial interval */
/*  f      function subprogram which evaluates f(x) for any x in */
/*         the interval  ax,bx */
/*  tol    desired length of the interval of uncertainty of the */
/*         final result ( .ge. 0.0d0) */


/*  output.. */

/*  zeroin abcissa approximating a zero of  f  in the interval ax,bx */


/*      it is assumed  that   f(ax)   and   f(bx)   have  opposite  signs */
/*  without  a  check.  zeroin  returns a zero  x  in the given interval */
/*  ax,bx  to within a tolerance  4*macheps*abs(x) + tol, where macheps */
/*  is the relative machine precision. */
/*      this function subprogram is a slightly  modified  translation  of */
/*  the algol 60 procedure  zero  given in  richard brent, algorithms for */
/*  minimization without derivatives, prentice - hall, inc. (1973). */



/*  compute eps, the relative machine precision */

    eps = 1.;
L10:
    eps /= 2.;
    tol1 = eps + 1.;
    if (tol1 > 1.) {
	goto L10;
    }

/* initialization */

    a = *ax;
    b = *bx;
    fa = (*f)(&a,&flag);
    if (flag) return 1.0;
    fb = (*f)(&b,&flag);
    if (flag) return 1.0;

/* begin step */

L20:
    c__ = a;
    fc = fa;
    d__ = b - a;
    e = d__;
L30:
    if (abs(fc) >= abs(fb)) {
	goto L40;
    }
    a = b;
    b = c__;
    c__ = a;
    fa = fb;
    fb = fc;
    fc = fa;

/* convergence test */

L40:
    tol1 = eps * 2. * abs(b) + *tol * .5;
    xm = (c__ - b) * (float).5;
    if (abs(xm) <= tol1) {
	goto L90;
    }
    if (fb == 0.) {
	goto L90;
    }

/* is bisection necessary */

    if (abs(e) < tol1) {
	goto L70;
    }
    if (abs(fa) <= abs(fb)) {
	goto L70;
    }

/* is quadratic interpolation possible */

    if (a != c__) {
	goto L50;
    }

/* linear interpolation */

    s = fb / fa;
    p = xm * 2. * s;
    q = 1. - s;
    goto L60;

/* inverse quadratic interpolation */

L50:
    q = fa / fc;
    r__ = fb / fc;
    s = fb / fa;
    p = s * (xm * 2. * q * (q - r__) - (b - a) * (r__ - 1.));
    q = (q - 1.) * (r__ - 1.) * (s - 1.);

/* adjust signs */

L60:
    if (p > 0.) {
	q = -q;
    }
    p = abs(p);

/* is interpolation acceptable */

    if (p * 2. >= xm * 3. * q - (d__1 = tol1 * q, abs(d__1))) {
	goto L70;
    }
    if (p >= (d__1 = e * .5 * q, abs(d__1))) {
	goto L70;
    }
    e = d__;
    d__ = p / q;
    goto L80;

/* bisection */

L70:
    d__ = xm;
    e = d__;

/* complete step */

L80:
    a = b;
    fa = fb;
    if (abs(d__) > tol1) {
	b += d__;
    }
    if (abs(d__) <= tol1) {
	b += d_sign(&tol1, &xm);
    }
    fb = (*f)(&b,&flag);
    if (flag) return 1.0;
    if (fb * (fc / abs(fc)) > 0.) {
	goto L20;
    }
    goto L30;

/* done */

L90:
    ret_val = b;
    return ret_val;
} /* zeroin_ */