/* @(#)Copyright (c), 1987, 1993 StatSci, Inc.  All rights reserved. */
static char whatssi[] = "@(#)zero_find.c version 3.14 created 3/12/93 ";

#include "S.h"

#define DEPS 1e-17
static char *func;
static double zfun(), zero();

zero_find(ff, guesses,tol)
char **ff; double *guesses, *tol;
{
	double t, fa,fb,a,b;
	func = *ff;
	t = *tol > 0. ? *tol : 100*DEPS;
	fa = zfun(a=guesses[0]);
	fb = zfun(b=guesses[1]);
	*guesses = t = zero(zfun,a,b,fa,fb,t);
	fprintf(stderr,"Answer: %lf\n",t);
}

static double zfun(z)
double z;
{
	char *mode = "double"; long length = 1;
	char *args[1]; double zz = z;
	args[0] = (char *)(&zz);
	call_S(func, 1L, args, &mode, &length, 0, 1, args);
	fprintf(stderr,"z: %lf, f: %lf\n",z,*((double *)args[0]));
	return(*((double *)args[0]));
}

static double zero(f,a,b,fa,fb,t)
/*  finds the real root of the function f lying between a and b */
/*  to within a tolerance of */
/*         6*PRECISION * abs(zero) + 2 * t */
/*  fa and fb must have opposite signs */
/*  this is brents algorithm */
/*  a, stored in sa, is the previous best approximation (i.e. the old b) */
/*  b, stored in sb, is the current best approximation */
/*  c is the most recently computed point satisfying (*f) (b)*f(c) .lt. 0 */
/*  d contains the correction to the approximation */
/*  e contains the previous value of d */
/*  m contains the bisection quantity (c-b)/2 */
double a,b,t, (*f)(), fa, fb;
{
#define abs(x) (x>0. ? x : -x)
	double tt,sa,sb,c,d,e,fc,tol,m,p,q,r,s;
	tt = t;
	sa = a;
	sb = b;
	if (fa==0.0) return(sa);
	if (fb==0.0) return(sb);
	if(fa*fb >0.0)
		Recover("zero - (*f) (a) and (*f) (b) are not of opposite sign",0);
	for(;;) {
		c = sa;
		fc = fa;
		e = sb-sa;
		d = e;
		do {
	/*  interchange b and c if abs (*f) (c) .lt. abs (*f) (b) */
			if (abs(fc)<abs(fb)) {
				sa = sb;
				sb = c;
				c = sa;
				fa = fb;
				fb = fc;
				fc = fa;
				}
			tol = 2.0*DEPS*abs(sb)+tt;
			m = 0.5*(c-sb);
	/*  success indicated by m reduced to under tolerance or */
	/*  by (*f) (b) = 0 */
			if (abs(m)<=tol||fb==0.0)return(sb);
	/*  a bisection is forced if e, the next-to-last correction */
	/*  was less than the tolerance or if the previous b gave */
	/*  a smaller (*f) (b).  otherwise go to 40. */
			if (abs(e)<tol||abs(fa)<abs(fb)) {
				e = m;
				d = e;
				}
			else {
				s = fb/fa;
	/*  quadratic interpolation can only be done if a (in sa) */
	/*  and c are different points. */
	/*  otherwise do the following linear interpolation */
				if (sa==c) {
					p = 2.0*m*s;
					q = 1.0-s;
					}
				else {
	/*  inverse quadratic interpolation */
					q = fa/fc;
					r = fb/fc;
					p = s*(2.0*m*q*(q-r)-(sb-sa)*(r-1.0));
					q = (q-1.0)*(r-1.0)*(s-1.0);
					}
				if (p<=0.0) p = -p;
				else q = -q;
	/*  update the quantities using the newly computed */
	/*  interpolate unless it would either force the */
	/*  new point too far to one side of the interval */
	/*  or would represent a correction greater than */
	/*  half the previous correction. */
	/*  in these last two cases - do the bisection */
				s = e;
				e = d;
				if (2.0*p<3.0*m*q-abs(tol*q)&&p<abs(0.5*s*q))
					d = p/q;
				else {
					e = m;
					d = e;
					}
				}
	/*  set a to the previous b */
			sa = sb;
			fa = fb;
	/*  if the correction to be made is smaller than */
	/*  the tolerance, just take a  delta step  (delta=tolerance) */
	/*         b = b + delta * sign(m) */
			if (abs(d)>tol) sb = sb+d;
			else if (m<=0.0) sb = sb-tol;
			else sb = sb+tol;
			fb = (*f) (sb);
	/*  if (*f) (b) and (*f) (c) have the same sign only */
	/*  linear interpolation (not inverse quadratic) */
	/*  can be done */
			}
			while( fc * fb < 0.0);
		}
}