in commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/GeneralizedContinuedFraction.java [364:449]
static double evaluate(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) {
// Relative error epsilon should not be zero to prevent drift in the event
// that the update ratio never achieves 1.0.
// Epsilon is the relative change allowed from 1. Configure the absolute limits so
// convergence requires: low <= deltaN <= high
// low = 1 - eps
// high = 1 / (1 - eps)
// High is always further from 1 than low in absolute distance. Do not store high
// but store the maximum absolute deviation from 1 for convergence = high - 1.
// If this is achieved a second check is made against low.
double low;
double eps;
if (epsilon > MIN_EPSILON && epsilon <= MAX_EPSILON) {
low = 1 - epsilon;
eps = 1 / low - 1;
} else {
// Precomputed defaults. Used when epsilon <= MIN_EPSILON
low = DEFAULT_LOW;
eps = DEFAULT_EPS;
}
double hPrev = updateIfCloseToZero(b0);
// Notes from Thompson and Barnett:
//
// Fraction convergent: hn = An / Bn
// A(-1) = 1, A0 = b0, B(-1) = 0, B0 = 1
// Compute the ratios:
// Dn = B(n-1) / Bn = 1 / (an * D(n-1) + bn)
// Cn = An / A(n-1) = an / C(n-1) + bn
//
// Ratio of successive convergents:
// delta n = hn / h(n-1)
// = Cn / Dn
// Avoid divisors being zero (less than machine precision) by shifting them to e.g. 1e-50.
double dPrev = 0.0;
double cPrev = hPrev;
for (int n = maxIterations; n > 0; n--) {
final Coefficient c = gen.get();
final double a = c.getA();
final double b = c.getB();
double dN = updateIfCloseToZero(b + a * dPrev);
final double cN = updateIfCloseToZero(b + a / cPrev);
dN = 1 / dN;
final double deltaN = cN * dN;
final double hN = hPrev * deltaN;
// If the fraction is convergent then deltaN -> 1.
// Computation of deltaN = 0 or deltaN = big will result in zero or overflow.
// Directly check for overflow on hN (this ensures the result is finite).
if (!Double.isFinite(hN)) {
throw new FractionException("Continued fraction diverged to " + hN);
}
// Check for underflow on deltaN. This allows fractions to compute zero
// if this is the convergent limit.
// Note: deltaN is only zero if dN > 1e-50 / min_value, or 2.02e273.
// Since dN is the ratio of convergent denominators this magnitude of
// ratio is a presumed to be an error.
if (deltaN == 0) {
throw new FractionException("Ratio of successive convergents is zero");
}
// Update from Thompson and Barnett to use <= eps in place of < eps.
// eps = high - 1
// A second check is made to ensure:
// low <= deltaN <= high
if (Math.abs(deltaN - 1) <= eps && deltaN >= low) {
return hN;
}
dPrev = dN;
cPrev = cN;
hPrev = hN;
}
throw new FractionException("Maximum iterations (%d) exceeded", maxIterations);
}