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;
}
}
}