""" brent(f, xl, xh, iterations = 100, ɛ = 1.0e-7) `xl` and `xh` bracket the root of f. """ function brent(f, xl, xh, iterations = 100, ɛ = 1.0e-7) a, fa = xl, f(xl) b, fb = xh, f(xh) @assert fa * fb < 0 "The root must be bracketed by `xl` and `xh`" c, fc = a, fa # c is the previous position d = e = b - c # d is the next step (b = b + d), e the previous one for i = 1:iterations if fb == 0 break # we have a root end if sign(fa) == sign(fb) # [a,b] should bracket the root a, fa = c, fc d = e = b - c end m = (a - b) / 2 # maximum error (half of the bracket interval) tol = 2 * ɛ * max(abs(b), 1) if abs(m) <= tol break # we are precise enough end if abs(e) >= tol && abs(fc) > abs(fb) # Use open methods local p, q # for efficient computation, d = p/q, p > 0 s = fb / fc if a == c # Secant method - when we have only 2 different points p = 2 * m * s q = 1 - s else # Inverse quadratic interpolation q, r = fc / fa, fb / fa p = s * (2 * m * q * (q - r) - (b - c) * (r - 1)) q = (q - 1) * (r - 1) * (s - 1) end if p > 0 q = -q else p = -p end if 2 * p < 3 * m * q - abs(tol * q) && p < abs(e * q / 2) d, e = p / q, d else d, e = m, m # Fall back to bisection when out of bounds end else # Use bracketing method [bisection] when bounds decrease too slowly d, e = m, m end c, fc = b, fb if abs(d) > tol b += d else b -= sign(b-a) * tol # Forces a minimal step end fb = f(b) end b end