private double brent()

in commons-numbers-rootfinder/src/main/java/org/apache/commons/numbers/rootfinder/BrentSolver.java [158:249]


    private double brent(DoubleUnaryOperator func,
                         double lo, double hi,
                         double fLo, double fHi) {
        double a = lo;
        double fa = fLo;
        double b = hi;
        double fb = fHi;
        double c = a;
        double fc = fa;
        double d = b - a;
        double e = d;

        final double t = absoluteAccuracy;
        final double eps = relativeAccuracy;

        while (true) {
            if (Math.abs(fc) < Math.abs(fb)) {
                a = b;
                b = c;
                c = a;
                fa = fb;
                fb = fc;
                fc = fa;
            }

            final double tol = 2 * eps * Math.abs(b) + t;
            final double m = 0.5 * (c - b);

            if (Math.abs(m) <= tol ||
                equalsZero(fb))  {
                return b;
            }
            if (Math.abs(e) < tol ||
                Math.abs(fa) <= Math.abs(fb)) {
                // Force bisection.
                d = m;
                e = d;
            } else {
                final double s = fb / fa;
                double p;
                double q;
                // The equality test (a == c) is intentional,
                // it is part of the original Brent's method and
                // it should NOT be replaced by proximity test.
                if (a == c) {
                    // Linear interpolation.
                    p = 2 * m * s;
                    q = 1 - s;
                } else {
                    // Inverse quadratic interpolation.
                    q = fa / fc;
                    final double r = fb / fc;
                    p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
                    q = (q - 1) * (r - 1) * (s - 1);
                }
                if (p > 0) {
                    q = -q;
                } else {
                    p = -p;
                }
                if (p >= 1.5 * m * q - Math.abs(tol * q) ||
                    p >= Math.abs(0.5 * e * q)) {
                    // Inverse quadratic interpolation gives a value
                    // in the wrong direction, or progress is slow.
                    // Fall back to bisection.
                    d = m;
                    e = d;
                } else {
                    e = d;
                    d = p / q;
                }
            }
            a = b;
            fa = fb;

            if (Math.abs(d) > tol) {
                b += d;
            } else if (m > 0) {
                b += tol;
            } else {
                b -= tol;
            }
            fb = func.applyAsDouble(b);
            if ((fb > 0 && fc > 0) ||
                (fb <= 0 && fc <= 0)) {
                c = a;
                fc = fa;
                d = b - a;
                e = d;
            }
        }
    }