static double evaluate()

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