static bool brent()

in include/maths/common/CSolvers.h [508:608]
85 lines of code
26 McCabe index (conditional complexity)

    static bool brent(double& a,
                      double& b,
                      double fa,
                      double fb,
                      const F& f,
                      std::size_t& maxIterations,
                      const EQUAL& equal,
                      double& bestGuess) {
        std::size_t n = std::max(maxIterations, std::size_t{1});
        if (fa == 0.0) {
            // Root at left bracket.
            bestGuess = b = a;
            maxIterations -= n;
            return true;
        }
        if (fb == 0.0) {
            // Root at right bracket.
            bestGuess = a = b;
            maxIterations -= n;
            return true;
        }
        if (std::signbit(fa) == std::signbit(fb)) {
            // Not bracketed.
            maxIterations -= n;
            return false;
        }

        if (std::fabs(fa) < std::fabs(fb)) {
            std::swap(a, b);
            std::swap(fa, fb);
        }

        bool bisected = true;
        double c = a;
        double fc = fa;
        double d = std::numeric_limits<double>::max();

        do {
            double s;
            bool canUseInverseQuadratic = ((fa != fc) && (fb != fc));
            bool canUseSecant = (fa != fb);
            if (canUseInverseQuadratic || canUseSecant) {
                s = canUseInverseQuadratic
                        ? inverseQuadraticInterpolate(a, b, c, fa, fb, fc)
                        : secantInterpolate(a, b, fa, fb);

                double e = (3.0 * a + b) / 4.0;

                if ((!(((s > e) && (s < b)) || ((s < e) && (s > b)))) ||
                    (bisected &&
                     ((std::fabs(s - b) >= std::fabs(b - c) / 2.0) || equal(b, c))) ||
                    (!bisected &&
                     ((std::fabs(s - b) >= std::fabs(c - d) / 2.0) || equal(c, d)))) {
                    // Use bisection.
                    s = bisect(a, b);
                    bisected = true;
                } else {
                    bisected = false;
                }
            } else {
                s = bisect(a, b);
                bisected = true;
            }

            double fs = f(s);
            shift(d, c, b);
            fc = fb;

            if (fs == 0.0) {
                // Root at s.
                bestGuess = a = b = s;
                maxIterations -= (n - 1);
                return true;
            }

            if (std::signbit(fa) == std::signbit(fs)) {
                a = s;
                fa = fs;
            } else {
                b = s;
                fb = fs;
            }

            if (std::fabs(fa) < std::fabs(fb)) {
                std::swap(a, b);
                std::swap(fa, fb);
            }
        } while (--n > 0 && !equal(a, b));

        if (b < a) {
            std::swap(a, b);
            std::swap(fa, fb);
        }
        bestGuess = (fa != fc) && (fb != fc)
                        ? inverseQuadraticInterpolate(a, b, c, fa, fb, fc)
                        : (fa != fb ? secantInterpolate(a, b, fa, fb) : bisect(a, b));
        bestGuess = std::min(std::max(a, bestGuess), b);
        maxIterations -= n;

        return true;
    }