Brent’s Root Finder
March 28, 2014
[ Today’s exercise was written by guest author Paul Hofstra, who holds a PhD in Theoretical Nuclear Physics from the Free University in Amsterdam and is now retired from a major oil company. Guest authors are always welcome; contact me at the link in the menu bar if you are interested in writing for Programming Praxis. ]
We discussed various algorithms for finding the roots of a function in two previous exercises. In the first of those exercises, we looked at the bisection method, which is “slow and steady” as it always finds a solution, and the secant method, which is sometimes very fast but also sometimes very slow. In the second of those two exercises, we looked at an algorithm of Theodorus Dekker that combines the bisection and secant methods in a way that retains the speed of the secant method where possible but retains the steadiness of the bisection method where necessary. Later, Richard Brent (the same Brent that improved Pollard’s rho algorithm for finding factors), observed that there are still some functions for which Dekker’s algorithm performs poorly, and he developed the method that is most commonly used today.
Brent augmented the linear interpolation of the secant method with a quadratic interpolation that tends to get closer to the root more quickly than the secant method, and he made conditions that force a bisection step at least every three steps, so his method can never be worse than three times the bisection method, and usually much better. Brent calculates a new point using both the secant and bisection methods, as does Dekker, but instead of simply accepting the secant method when it is better, he calculates another potential point using the quadratic method and takes that instead if it is better than the secant point. You can read Brent’s description of his algorithm at http://maths-people.anu.edu.au/~brent/pd/rpb005.pdf.
Your task is to implement Brent’s root finder. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.