roots.cpp

Go to the documentation of this file.

/// @file roots.cpp
#include "analysis/roots.hpp"
#include <cmath>
#include <stdexcept>
 
namespace num {
 
RootResult bisection(ScalarFn f, real a, real b, real tol, idx max_iter) {
    real fa = f(a), fb = f(b);
    if (fa * fb > 0.0)
        throw std::invalid_argument("bisection: f(a) and f(b) must have opposite signs");
 
    for (idx i = 0; i < max_iter; ++i) {
        real mid = 0.5 * (a + b);
        real fm  = f(mid);
        if (std::abs(fm) < tol || 0.5 * (b - a) < tol)
            return {mid, i + 1, std::abs(fm), true};
        if (fa * fm < 0.0) { b = mid; fb = fm; }
        else               { a = mid; fa = fm; }
    }
    real mid = 0.5 * (a + b);
    return {mid, max_iter, std::abs(f(mid)), false};
}
 
RootResult newton(ScalarFn f, ScalarFn df, real x0, real tol, idx max_iter) {
    real x = x0;
    for (idx i = 0; i < max_iter; ++i) {
        real fx  = f(x);
        if (std::abs(fx) < tol)
            return {x, i + 1, std::abs(fx), true};
        real dfx = df(x);
        if (std::abs(dfx) < 1e-14)
            return {x, i + 1, std::abs(fx), false};  // near-zero derivative
        x -= fx / dfx;
    }
    return {x, max_iter, std::abs(f(x)), false};
}
 
RootResult secant(ScalarFn f, real x0, real x1, real tol, idx max_iter) {
    real f0 = f(x0), f1 = f(x1);
    for (idx i = 0; i < max_iter; ++i) {
        if (std::abs(f1) < tol)
            return {x1, i + 1, std::abs(f1), true};
        real df = f1 - f0;
        if (std::abs(df) < 1e-14)
            return {x1, i + 1, std::abs(f1), false};  // stagnation
        real x2 = x1 - f1 * (x1 - x0) / df;
        x0 = x1; f0 = f1;
        x1 = x2; f1 = f(x1);
    }
    return {x1, max_iter, std::abs(f1), false};
}
 
RootResult brent(ScalarFn f, real a, real b, real tol, idx max_iter) {
    real fa = f(a), fb = f(b);
    if (fa * fb > 0.0)
        throw std::invalid_argument("brent: f(a) and f(b) must have opposite signs");
 
    // c is the point such that [b,c] is always a valid bracket
    real c = a, fc = fa;
    real d = b - a, e = d;
 
    for (idx i = 0; i < max_iter; ++i) {
        // ensure b is the best estimate
        if (fb * fc > 0.0) { c = a; fc = fa; d = e = b - a; }
        if (std::abs(fc) < std::abs(fb)) {
            a = b; fa = fb;
            b = c; fb = fc;
            c = a; fc = fa;
        }
 
        real tol1 = 2.0 * 1e-15 * std::abs(b) + 0.5 * tol;
        real mid  = 0.5 * (c - b);
 
        if (std::abs(mid) <= tol1 || std::abs(fb) < tol)
            return {b, i + 1, std::abs(fb), true};
 
        if (std::abs(e) >= tol1 && std::abs(fa) > std::abs(fb)) {
            // attempt interpolation
            real s = fb / fa;
            real p, q;
            if (a == c) {
                // secant step
                p = 2.0 * mid * s;
                q = 1.0 - s;
            } else {
                // inverse quadratic interpolation
                real r = fb / fc;
                real t = fa / fc;
                p = s * (2.0 * mid * t * (t - r) - (b - a) * (r - 1.0));
                q = (t - 1.0) * (r - 1.0) * (s - 1.0);
            }
            if (p > 0.0) q = -q; else p = -p;
 
            // accept interpolation only if step is smaller than previous
            real e_prev = e;
            if (2.0 * p < std::min(3.0 * mid * q - std::abs(tol1 * q),
                                   std::abs(e_prev * q))) {
                e = d;
                d = p / q;
            } else {
                d = mid; e = mid;
            }
        } else {
            d = mid; e = mid;
        }
 
        a = b; fa = fb;
        b += (std::abs(d) > tol1) ? d : (mid > 0.0 ? tol1 : -tol1);
        fb = f(b);
    }
    return {b, max_iter, std::abs(fb), false};
}
 
} // namespace num