Brent's Zero Finding Algorithm  (also called ZEROIN)
(A search for Brent or Zeroin should find other versions of this algorithm)

Given a, b  with f(a)f(b) < 0, tolerance delta and unit roundoff eps
find a root x within 2 delta of real root of f

init:	fb = f(b)	% b is the best approximation to the root
	fa = f(a)	% a is the last best approximation
setup:	c = a		% the root is always between b and c
	fc = fa
	e = b - a	% e is the last correction
	d = e		% d is the latest correction
order:	if |fc|<|fb| then	% always must have |fb| < |fc|
	   a = b, b = c, c = a
	   fa = fb, fb = fc, fc = fa
	endif
main:	t = 2 eps |b| + delta	% actual tolerance, accounts for roundoff
	m = (c-b)/2		% b + m is bisection approximation
	if (|m|<t) or (fb=0) goto end:	% within t of root
	if (|e|<t) or (|fa|<|fb|) then	% take bisection
	   e = m
	   d = e
	else				% try interpolation
	   s = fb/fa			% the interpolated value is
	   if (a=c) then		%    i = b + p/q
		p = 2ms	
		q = 1-s
	   else
		q = fa/fc
		r = fb/fc
		p = s(2mq(q-r) - (b-a)(r-1))
		q = (q-1)(r-1)(s-1)
	   endif
	   if (p>0) then	% always make p>0
		q = -q
	   else
		p = -p
	   endif
	   s = e
	   e = d		% save last correction
%  we only use the interpolated value if it (1) is between b and c, and
%                                           (2) we are converging fast enough
	   if (2p<3mq-|tq|) and (p<|sq/2|) then 	% it is good
		d = p/q
	   else						% bisection is better
		e = m
		d = e
	   endif
update:	   a = b
	   fa = fb
	   if (|d|>t) then	% good step
		b = b + d
	   else			% take a step of size t
		if (m>0) then
		   b = b + t
		else
		   b = b - t
		endif
	   endif
	   fb = f(b)
	   if (fb fc >0) goto setup:
	   goto order:
	endif
end:	x = b