Hurmet functions for numerically finding roots to equations

rootFinding.txt

# Three Hurmet functions for numerically finding roots of equations.

# First, bisection. The first argument is the function, f.
# Select a and b such that f(a) and f(b) have opposite signs.
# The optional argument ε enables you to define a desired precision.

function bisection(f, a, b; ε = 1e-15)
    fa = f(a)
    fb = f(b)
    if fa × fb > 0 throw "Error. Invalid starting bracket."
    while true
        x = (a + b) / 2
        fx = f(x)
        if |fx| ≤ ε return x
        if sign(fa) == sign(fx)
            a = x
        else
            b = x
        end
    end
end

# Next, Newton's method. Fast when it works, but sometimes does not converge.
# Don’t use it on a periodic function. It will lock up a browser tab.
# Input the function f, its first derivative f′, and a starting guess. The ε is optional.

function newton(f, fPrime, guess; ε = 1e-15)
    x = guess
    while true
        fx = f(x)
        if |fx| ≤ ε return x
        x = x - fx / fPrime(x)
    end
end

# Finally, Brent's method. Faster than bisection. Sure to find a result, so long
# as the function is continuous.
# Select a and b such that f(a) and f(b) have opposite signs.
# Adapted from John D Cook: Three Methods for Root-finding in C#
#    https://www.codeproject.com/Articles/79541/Three-Methods-for-Root-finding-in-C

function brent(f, a, b; ε = 1e-15)
    fa = f(a)
    fb = f(b)
    if fa × fb > 0 throw "Error. Invalid starting bracket."
    c = a
    fc = fa
    e = b - a
    d = e
    while true
        if |fb| ≤ ε return b
        if (fb > 0.0 and fc > 0.0) or (fb ≤ 0.0 and fc ≤ 0.0)
            c = a
            fc = fa
            e = b - a
            d = e
        end
        if |fc| < |fb|
            a = b
            b = c
            c = a
            fa = fb
            fb = fc
            fc = fa
        end
        tol = 2 ε · |b| + ε
        m = 0.5 · (c - b)  # error estimate
        if |m| > tol and fb ≠ 0.0
            if |e| < tol or |fa| ≤ |fb|
                # Use bisection
                e = m
                d = e
            else
                s = fb / fa
                if a == c
                    # linear intepolation
                    p = 2 m s
                    q = 1.0 - s
                else
                    # Inverse quadratic interpolation
                    q = fa / fc
                    r = fb / fc
                    p = s * (2 m q * (q - r) - (b - a) * (r - 1.0))
                    q = (q - 1.0) * (r - 1.0) * (s - 1.0)
                end
                if p > 0.0
                    q = -q
                else
                    p = -p
                end
                s = e
                e = d
                if 2 p < 3 m q - |tol * q| and p < |0.5 s q|
                    d = p / q
                else
                    e = m
                    d = e
                end
            end
            a = b
            fa = fb
            if |d| > tol
                b = b + d
            elseif m > 0.0
                b = b + tol
            else
                b = b - tol
            end
            fb = f(b)
        else
          return b
        end
    end
end