/*
* Brent's method. Returns false if bad params. Returns [x,y] of solution
* otherwise.
* Taken from Numerical Recipes by Press, et al.
*/
brent(brent_func_name,x1,x2,tol):=
block([a:x1,b:x2,c,d,e,eps:1.b1^(-bfprecision+1),
fa,fb,fc,iter,p,q,r,s,tol1,xm],
fa:brent_func_name(a),
fb:brent_func_name(b),
num_brent_func_evals:2,
if fb*fa > 0.b0 then return(false),
fc:fb,
return (for iter:1 step 1 thru brent_itmax do (
if (bent_brent_verbose) then (
print("a", iter, a, fa, b, fb, c, fc),
0),
if fb*fc > 0.b0 then ( c:a, fc:fa, d:b-a, e:d ),
if (bent_brent_verbose) then (
print("b", iter, a, fa, b, fb, c, fc),
0),
if abs(fc) < abs(fb) then (
a:b, b:c, c:a, fa:fb, fb:fc, fc:fa,
0),
if (bent_brent_verbose) then (
print("c", iter, a, fa, b, fb, c, fc),
0),
tol1:2.b0*eps*abs(b)+0.5b0*tol,
xm:0.5b0*(c-b),
if (bent_brent_verbose) then (
print("d", iter, "tol1 = ", tol1, "xm = ", xm),
0),
if abs(xm) <= tol1 or zerop(fb) then return([b,fb]),
if abs(e) >= tol1 and abs(fa) > abs(fb) then (
s:fb/fa,
if a = c then (
p:2.b0*xm*s,
q:1.b0-s,
if (bent_brent_verbose) then (
print("e", iter, "s = ", s,
"p = ", p, "q = ", q),
0),
0) else (
q:fa/fc,
r:fb/fc,
p:s*(2.b0*xm*q*(q-r)-(b-a)*(r-1.b0)),
q:(q-1.b0)*(r-1.b0)*(s-1.b0),
if (bent_brent_verbose) then (
print("f", iter, "s = ", s,
"p = ", p, "q = ", q,
"r = ", r),
0),
1),
if p > 0.b0 then (
q:-q,
0),
p:abs(p),
if 2.b0*p < min(3.b0*xm*q-abs(tol1*q),abs(e*q)) then (
e:d,
d:p/q,
if (bent_brent_verbose) then (
print("t", iter, "interp",
"d = ", d, "e = ", e),
0),
2) else (
d:xm,
e:d,
if (bent_brent_verbose) then (
print("u", iter, "interp failed",
"d = ", d, "e = ", e),
0),
0),
3) else (
d:xm,
e:d,
if (bent_brent_verbose) then (
print("v", iter, "interp too slow",
"d = ", d, "e = ", e),
0),
0),
a:b,
fa:fb,
if abs(d) > tol1 then (
b:b+d,
if (bent_brent_verbose) then (
print("w", iter, "accept interp/bisect",
"adj = ", d),
0),
0) else if xm > 0 then (
b:b+abs(tol1),
if (bent_brent_verbose) then (
print("x", iter, "reject interp/bisect",
"adj = ", abs(tol1)),
0),
0) else (
b:b-abs(tol1),
if (bent_brent_verbose) then (
print("y", iter, "reject interp/bisect",
"adj = ", -abs(tol1)),
0),
0),
fb:brent_func_name(b),
num_brent_func_evals:num_brent_func_evals + 1,
4)),
5)$
brent_itmax:100$