double x0 = min;
double x1 = max;
double y0 = f.value(x0);
double y1 = f.value(x1);
if ((y0 > 0) == (y1 > 0)) {
throw new MathException("Interval doesn't bracket a zero.");
}
double x2 = x0;
double y2 = y0;
double delta = x1 - x0;
double oldDelta = delta;
int i = 0;
while (i < maximalIterationCount) {
if (Math.abs(y2) < Math.abs(y1)) {
x0 = x1;
x1 = x2;
x2 = x0;
y0 = y1;
y1 = y2;
y2 = y0;
}
if (Math.abs(y1) <= functionValueAccuracy) {
// Avoid division by very small values. Assume
// the iteration has converged (the problem may
// still be ill conditioned)
setResult(x1, i);
return result;
}
double dx = (x2 - x1);
double tolerance =
Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy);
if (Math.abs(dx) <= tolerance) {
setResult(x1, i);
return result;
}
if ((Math.abs(oldDelta) < tolerance) ||
(Math.abs(y0) <= Math.abs(y1))) {
// Force bisection.
delta = 0.5 * dx;
oldDelta = delta;
} else {
double r3 = y1 / y0;
double p;
double p1;
if (x0 == x2) {
// Linear interpolation.
p = dx * r3;
p1 = 1.0 - r3;
} else {
// Inverse quadratic interpolation.
double r1 = y0 / y2;
double r2 = y1 / y2;
p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
}
if (p > 0.0) {
p1 = -p1;
} else {
p = -p;
}
if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) ||
p >= Math.abs(0.5 * oldDelta * p1)) {
// Inverse quadratic interpolation gives a value
// in the wrong direction, or progress is slow.
// Fall back to bisection.
delta = 0.5 * dx;
oldDelta = delta;
} else {
oldDelta = delta;
delta = p / p1;
}
}
// Save old X1, Y1
x0 = x1;
y0 = y1;
// Compute new X1, Y1
if (Math.abs(delta) > tolerance) {
x1 = x1 + delta;
} else if (dx > 0.0) {
x1 = x1 + 0.5 * tolerance;
} else if (dx <= 0.0) {
x1 = x1 - 0.5 * tolerance;
}
y1 = f.value(x1);
if ((y1 > 0) == (y2 > 0)) {
x2 = x0;
y2 = y0;
delta = x1 - x0;
oldDelta = delta;
}
i++;
}
throw new MathException("Maximal iteration number exceeded.");
}