include	<mach.h>
define	ITMAX		100
define	FACTOR		1.6
define	NTRY		50

# BRENT -- Find the roots of a function by Brent's method

# This procedure uses Brent's method to find one root of a function. Brent's
# method uses inverse quadratic interpolation suplemented by bisection when
# the trial root lies outside the bracket. The code was converted from the
# Fortran routine ZBRENT in Numerical Recipes, section 9.3.
#
# B.Simon	11-Sep-90	First Code

real procedure brent (func, x1, x2, tol)

real	func		# i: Function for which the root is sought
real	x1		# i: First point in bracket around root
real	x2		# i: Second point in bracket around root (x1 != x2)
real	tol		# i: Accuracy required in root
#--
extern	func

int	iter
real	a, b, c, d, e, fa, fb, fc
real	xm, p, q, r, s, tol1

bool	bracket_root()

errchk	func, bracket_root

begin
	a = x1
	b = x2
	fa = func (a)
	fb = func (b)

	# Check to see if the bracket is valid
	# If not, search for a new bracket

	if (fb * fa > 0.) {
	    if (! bracket_root (func, a, b)) {
		call error (1, "Root must be bracketed for brent method")
	    } else {
		fa = func (a)
		fb = func (b)
	    }
	}

	fc = fb
	do iter = 1, ITMAX {

	    # Rename A,B,C and adjust bounding interval D

	    if (fb * fc > 0.) {
		c = a
		fc = fa
		d = b - a
		e = d
	    }

	    if(abs(fc) < abs(fb)) {
		a = b
		b = c
		c = a
		fa = fb
		fb = fc
		fc = fa
	    }

	    # Convergence check and successful return from function

	    tol1 = 2. * EPSILONR * abs(b) + 0.5 * tol
	    xm = .5 * (c - b)
	    if(abs(xm) <= tol1 || fb == 0.)
		return (b)

	    if(abs(e) >= tol1 && abs(fa) > abs(fb)) {
		# Attempt inverse quadratic interpolation

		s = fb / fa
		if(a == c) {
		    p = 2. * xm * s
		    q = 1. - s
		} else {
		    q = fa / fc
		    r = fb / fc
		    p = s * (2. * xm * q * (q - r) - (b - a) * (r - 1.))
		    q = (q - 1.) * (r - 1.) * (s - 1.)
		}

		# Check whether in bounds

		if (p > 0.) 
		    q = -q

		p = abs(p)
		if(2. * p < min(3. * xm * q - abs(tol1 * q), abs(e * q))) {
		    # Accept interpolation

		    e = d
		    d = p / q

		} else {
		    # Interpolation failed, use bisection

		    d = xm
		    e = d
		}

	    } else {
		# Bounds decreasing too slowly, use bisection

		d = xm
		e = d
	    }

	    # Move last best guess to A

	    a = b
	    fa = fb

	    # Evaluate new trial root

	    if(abs(d) > tol1) {
		b = b + d
	    } else {
		b = b + sign (tol1, xm)
	    }
	    fb = func (b)
	}

	# Error exit, no convergence to root

	call error (1, "Brent exceeding maximum iterations")
	return (b)
end

# BRACKET_ROOT -- Find an initial bracket for a root

# This procedure takes an initial guessed bracket and expands the range
# until it contains a root. Since the method is not guaranteed to find
# a bracket, it returns true if the bracket is found and false otherwise.
# The function is taken from the Fortran routine ZBRAC in Numerical 
# Recipes, section 9.1

bool procedure bracket_root (func, x1, x2)

real	func		#  i: Function for which bracket is sought
real	x1		# io: First point in bracket around root
real	x2		# io: Second point in bracket around root
#--
extern	func

bool	success
int	j
real	f1, f2

begin

	if (x1 == x2)
	    call error (1, "bracket_root: You have to guess an initial range")

	f1 = func(x1)
	f2 = func(x2)
	success = false

	do j = 1, NTRY {
	    if (f1 * f2 < 0.) {
		success = true
		break
	    }

	    if (abs(f1) < abs(f2)) {
		x1 = x1 + FACTOR * (x1 - x2)
		f1 = func (x1)
	    } else {
		x2 = x2 + FACTOR * (x2 - x1)
		f2 = func (x2)
	    }
	}

	return (success)

end