Brent'ss Algorithm
Introduction
Brent’s algorithm, also known as Brent’s method or the Brent–dekker method, is a rootfinding algorithm that combines the bisection method, the secant method, and inverse quadratic interpolation. It is an iterative algorithm that can find the root of a continuous function in a given interval. The algorithm is widely used in numerical analysis, optimization, and computer science.
Implementation
The following is an implementation of Brent’s algorithm in OCaml. The function takes a function f
and two endpoints a
and b
of an interval [a, b]
as input, and returns an approximate root of f
in the interval. The optional parameters max_iter
and xtol
are the maximum number of iterations and the tolerance for convergence, respectively.
let brent ?(max_iter = 1000) ?(xtol = 1e6) f a b =
let fa = f a in
let fb = f b in
let error () =
let s = Printf.sprintf "f(a) *. f(b) = %g *. %g should be negative." fa fb in
Owl_exception.INVALID_ARGUMENT s
in
Owl_exception.verify (fa *. fb < 0.) error;
if fa = 0.
then a
else if fb = 0.
then b
else (
let xa = ref a in
let xb = ref b in
let xc = ref b in
let fc = ref fb in
let fa = ref fa in
let fb = ref fb in
let d = ref infinity in
let e = ref infinity in
let p = ref infinity in
let q = ref infinity in
let r = ref infinity in
let eps = 3e16 in
try
for _i = 1 to max_iter do
if (!fb > 0. && !fc > 0.)  (!fb < 0. && !fc < 0.)
then (
xc := !xa;
fc := !fa;
d := !xb . !xa;
e := !d);
if abs_float !fc < abs_float !fb
then (
xa := !xb;
xb := !xc;
xc := !xa;
fa := !fb;
fb := !fc;
fc := !fa);
let tol = (2. *. eps *. abs_float !xb) +. (0.5 *. xtol) in
let xm = 0.5 *. (!xc . !xb) in
assert (abs_float xm >= tol && !fb != 0.);
(* 1st strategy: inverse quadratic interpolation *)
if abs_float !e >= tol && abs_float !fa > abs_float !fb
then (
let s = !fb /. !fa in
if !xa = !xc
then (
p := 2. *. xm *. s;
q := 1. . s)
else (
q := !fa /. !fc;
r := !fb /. !fc;
p := s *. ((2. *. xm *. !q *. (!q . !r)) . ((!xb . !xa) *. (!r . 1.)));
q := (!q . 1.) *. (!r . 1.) *. (s . 1.));
if !p > 0. then q := . !q;
p := abs_float !p;
let min1 = (3. *. xm *. !q) . abs_float (tol *. !q) in
let min2 = abs_float (!e *. !q) in
if 2. *. !p < min min1 min2
then (
e := !d;
d := !p /. !q)
else (
d := xm;
e := !d)
(* 2nd strategy: bisection method *))
else (
d := xm;
e := !d);
(* adjust the position *)
xa := !xb;
fa := !fb;
if abs_float !d > tol
then xb := !xb +. !d
else xb := !xb +. if tol > 0. then xm else .xm;
fb := f !xb
done;
!xb
with
 _ > !xb)
Here is an example of using the function to find the root of the function f(x) = x^3  2x  5
in the interval [2, 3]
.
let f x = x ** 3. . 2. *. x . 5.
let root = brent f 2. 3.
let () = Printf.printf "The root of f(x) = x^3  2x  5 in [2, 3] is %g.\n" root
The output should be:
The root of f(x) = x^3  2x  5 in [2, 3] is 2.09455.
Stepbystep Explanation
Brent’s algorithm is an iterative algorithm that combines the bisection method, the secant method, and inverse quadratic interpolation. The algorithm starts with two endpoints a
and b
of an interval [a, b]
that contains a root of a continuous function f(x)
. The algorithm then iteratively refines the estimate of the root until it converges to a desired level of accuracy.
 Initialize variables
Let fa
and fb
be the values of f
at a
and b
, respectively. If fa
or fb
is zero, return the corresponding endpoint as the root. Otherwise, initialize the following variables:
xa
: the previous value ofxb
xb
: the current estimate of the rootxc
: the previous value ofxa
fa
: the value off
ata
fb
: the value off
atb
fc
: the value off
atc
d
: the step sizee
: the previous value ofd
p
,q
,r
: variables used in inverse quadratic interpolationeps
: a small number used for floatingpoint comparison
 Check for sign change
Verify that fa
and fb
have opposite signs. If not, raise an error.
 Iterate until convergence
Repeat the following steps until convergence or the maximum number of iterations is reached:

Check for convergence
Compute the tolerance
tol
for convergence. Iff(b)
is small enough or the step sized
is smaller thantol
, returnb
as the root. 
Compute the step size
Compute the midpoint
xm
of the interval[b, c]
, wherec
is the previous value ofa
orb
. Ife >= tol
andf(a) > f(b)
, use inverse quadratic interpolation to compute a new step sized
. Otherwise, use bisection method to computed
. 
Update the estimate of the root
Compute the new estimate of the root by adding
d
toxb
. Ifd
is too small, add or subtracttol
instead. 
Evaluate the function
Compute the value of
f
at the new estimatexb
. 
Update variables
Update the variables
xa
,xb
,xc
,fa
,fb
,fc
,d
, ande
with their new values. 
Check for convergence
If
f(b)
is small enough or the step sized
is smaller thantol
, returnb
as the root. 
Check for oscillation
If the function value at the new estimate
xb
has the same sign as the function value at the previous estimatexc
, replacexc
andfc
withxa
andfa
, and setd
ande
to their initial values. 
Check for convergence
If
f(b)
is small enough or the step sized
is smaller thantol
, returnb
as the root. 
Choose the interpolation method
If
f(b) < f(a)
, use secant method to compute a new step sized
. Otherwise, use inverse quadratic interpolation. 
Update the estimate of the root
Compute the new estimate of the root by adding
d
toxb
. Ifd
is too small, add or subtracttol
instead. 
Evaluate the function
Compute the value of
f
at the new estimatexb
. 
Update variables
Update the variables
xa
,xb
,xc
,fa
,fb
,fc
,d
, ande
with their new values. 
Return the root
If the algorithm does not converge within the maximum number of iterations, return the last estimate of the root.
Complexity Analysis
The time complexity of Brent’s algorithm is O(log(tol/eps) * f_evals), where tol is the tolerance for convergence, eps is a small number used for floatingpoint comparison, and f_evals is the number of evaluations of the function f
required by the algorithm. The algorithm is guaranteed to converge if the function f
is continuous and has a root in the initial interval [a, b]
. The algorithm is also robust to some types of singularities, such as poles and branch cuts. However, the algorithm may fail to converge if the function has multiple roots or if the initial interval is too large. In practice, Brent’s algorithm is often more efficient than other rootfinding algorithms, such as the bisection method and the secant method, especially for functions that are smooth but not necessarily analytic.