Search code examples
calgorithmimplementationfortrannumerical

What is "e" variable in popular implementations of Brent root finding algorithm?


I am reading the standard (Numerical Recipes and GSL C versions are identical) implementation of Brent root finding algorithm, and cannot understand the meaning of variable "e". The usage suggests that "e" is supposed to be the previous distance between the brackets. But then, why is it set to "xm" (half the distance) when we use bisection?


Solution

  • I'm not familiar with the algorithm. However, I can compare the C source and the Wikipedia description of the algorithm. The algorithm seems straight forward-ish (if you're familiar with methods to find roots), but the C implementation looks like a direct port of the fortran, so it's rather hard to read.

    My best guess is that e is related to the loop conditional.

    Wikipedia says (line 8 of the algorithm): repeat until f(b or s) = 0 or |b − a| is small enough (convergence)

    The C source says: e = b - a, then later if (fabs(e) <= tol ....

    I'd hope that the purpose of the variables would be described clearly in the book, but apparently not :)


    Ok, here you go. I found the original implementation (in algol 60) here. In addition to a nice description of the algorithm, it says (starting on page 50):

    let e be the value of p/q at the step before the last one. If |e|< δ or|p/q|1/2|e| then do a bisection, otherwise we do either a bisection or interpolation just as in Dekker's algorithm. Thus |e| decreases by at least a factor of two on every second step, and when |e|< δ a bisection must be done. (After a bisection we take e = m for the next step.)

    So the addition of e is Brent's "main modification" of Dekker's algorithm.