GeneralizedContinuedFraction.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.numbers.fraction;
import java.util.function.Supplier;
/**
* Provides a means to evaluate
* <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">generalized continued fractions</a>.
*
* <p>The continued fraction uses the following form for the numerator ({@code a}) and
* denominator ({@code b}) coefficients:
* <pre>
* a1
* b0 + ------------------
* b1 + a2
* -------------
* b2 + a3
* --------
* b3 + ...
* </pre>
*
* <p>A generator of the coefficients must be provided to evaluate the continued fraction.
*
* <p>The implementation of the fraction evaluation is based on the modified Lentz algorithm
* as described on page 508 in:
*
* <ul>
* <li>
* I. J. Thompson, A. R. Barnett (1986).
* "Coulomb and Bessel Functions of Complex Arguments and Order."
* Journal of Computational Physics 64, 490-509.
* <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
* https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
* </li>
* </ul>
*
* @see <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">Wikipedia: Generalized continued fraction</a>
* @see <a href="https://en.wikipedia.org/wiki/Generalized_continued_fraction">MathWorld: Generalized continued fraction</a>
* @since 1.1
*/
public final class GeneralizedContinuedFraction {
/**
* The value for any number close to zero.
*
* <p>"The parameter small should be some non-zero number less than typical values of
* eps * |b_n|, e.g., 1e-50".
*/
static final double SMALL = 1e-50;
/** Default maximum number of iterations. */
static final int DEFAULT_ITERATIONS = Integer.MAX_VALUE;
/**
* Minimum relative error epsilon. Equal to 1 - Math.nextDown(1.0), or 2^-53.
*
* <p>The epsilon is used to compare the change in the magnitude of the fraction
* convergent to 1.0. In theory eps can be 2^-53 reflecting the smallest reduction in
* magnitude possible i.e. {@code next = previous * Math.nextDown(1.0)}, or zero
* reflecting exact convergence.
*
* <p>If set to zero then the algorithm requires exact convergence which may not be possible
* due to floating point error in the algorithm. For example the golden ratio will not
* converge.
*
* <p>The minimum value will stop the recursive evaluation at the smallest possible
* increase or decrease in the convergent.
*/
private static final double MIN_EPSILON = 0x1.0p-53;
/** Maximum relative error epsilon. This is configured to prevent incorrect usage. Values
* higher than 1.0 invalidate the relative error lower bound of {@code (1 - eps) / 1}.
* Set to 0.5 which is a very weak relative error tolerance. */
private static final double MAX_EPSILON = 0.5;
/** Default low threshold for change in magnitude. Precomputed using MIN_EPSILON.
* Equal to 1 - 2^-53. */
private static final double DEFAULT_LOW = 1 - MIN_EPSILON;
/** Default absolute difference threshold for change in magnitude. Precomputed using MIN_EPSILON.
* Equal to {@code 1 / (1 - 2^-53) = 2^-52}. */
private static final double DEFAULT_EPS = 0x1.0p-52;
/**
* Defines the <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">
* {@code n}-th "a" and "b" coefficients</a> of the continued fraction.
*
* @since 1.1
*/
public static final class Coefficient {
/** "a" coefficient. */
private final double a;
/** "b" coefficient. */
private final double b;
/**
* @param a "a" coefficient
* @param b "b" coefficient
*/
private Coefficient(double a, double b) {
this.a = a;
this.b = b;
}
/**
* Returns the {@code n}-th "a" coefficient of the continued fraction.
*
* @return the coefficient <code>a<sub>n</sub></code>.
*/
public double getA() {
return a;
}
/**
* Returns the {@code n}-th "b" coefficient of the continued fraction.
*
* @return the coefficient <code>b<sub>n</sub></code>.
*/
public double getB() {
return b;
}
/**
* Create a new coefficient.
*
* @param a "a" coefficient
* @param b "b" coefficient
* @return the coefficient
*/
public static Coefficient of(double a, double b) {
return new Coefficient(a, b);
}
}
/** No instances. */
private GeneralizedContinuedFraction() {}
/**
* Evaluates the continued fraction.
*
* <p>Note: The first generated partial numerator a<sub>0</sub> is discarded.
*
* @param gen Generator of coefficients.
* @return the value of the continued fraction.
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number of
* iterations is reached before the expected convergence is achieved.
* @see #value(Supplier,double,int)
*/
public static double value(Supplier<Coefficient> gen) {
return value(gen, MIN_EPSILON, DEFAULT_ITERATIONS);
}
/**
* Evaluates the continued fraction.
*
* <p>Note: The first generated partial numerator a<sub>0</sub> is discarded.
*
* @param gen Generator of coefficients.
* @param epsilon Maximum relative error allowed.
* @return the value of the continued fraction.
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number of
* iterations is reached before the expected convergence is achieved.
* @see #value(Supplier,double,int)
*/
public static double value(Supplier<Coefficient> gen, double epsilon) {
return value(gen, epsilon, DEFAULT_ITERATIONS);
}
/**
* Evaluates the continued fraction.
* <pre>
* a1
* b0 + ------------------
* b1 + a2
* -------------
* b2 + a3
* --------
* b3 + ...
* </pre>
*
* <p>Setting coefficient a<sub>n</sub> to zero will signal the end of the recursive evaluation.
*
* <p>Note: The first generated partial numerator a<sub>0</sub> is discarded.
*
* <p><b>Usage Note</b>
*
* <p>This method is not functionally identical to calling
* {@link #value(double, Supplier, double, int)} with the generator configured to
* provide coefficients from n=1 and supplying b<sub>0</sub> separately. In some cases
* the computed result from the two variations may be different by more than the
* provided epsilon. The other method should be used if b<sub>0</sub> is zero or very
* small. See the corresponding javadoc for details.
*
* @param gen Generator of coefficients.
* @param epsilon Maximum relative error allowed.
* @param maxIterations Maximum number of iterations.
* @return the value of the continued fraction.
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number of
* iterations is reached before the expected convergence is achieved.
* @see #value(double, Supplier, double, int)
*/
public static double value(Supplier<Coefficient> gen, double epsilon, int maxIterations) {
// Use the first b coefficient to seed the evaluation of the fraction.
// Coefficient a is discarded.
final Coefficient c = gen.get();
return evaluate(c.getB(), gen, epsilon, maxIterations);
}
/**
* Evaluates the continued fraction.
*
* <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
* Both of the first generated terms a and b are used. This fraction evaluation
* can be used when:
* <ul>
* <li>b<sub>0</sub> is not part of a regular series
* <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component
* <li>b<sub>0</sub> is very small and the result is expected to approach zero
* </ul>
*
* @param b0 Coefficient b<sub>0</sub>.
* @param gen Generator of coefficients.
* @return the value of the continued fraction.
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number
* of iterations is reached before the expected convergence is achieved.
* @see #value(double,Supplier,double,int)
*/
public static double value(double b0, Supplier<Coefficient> gen) {
return value(b0, gen, MIN_EPSILON, DEFAULT_ITERATIONS);
}
/**
* Evaluates the continued fraction.
*
* <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
* Both of the first generated terms a and b are used. This fraction evaluation
* can be used when:
* <ul>
* <li>b<sub>0</sub> is not part of a regular series
* <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component
* <li>b<sub>0</sub> is very small and the result is expected to approach zero
* </ul>
*
* @param b0 Coefficient b<sub>0</sub>.
* @param gen Generator of coefficients.
* @param epsilon Maximum relative error allowed.
* @return the value of the continued fraction.
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number
* of iterations is reached before the expected convergence is achieved.
* @see #value(double,Supplier,double,int)
*/
public static double value(double b0, Supplier<Coefficient> gen, double epsilon) {
return value(b0, gen, epsilon, DEFAULT_ITERATIONS);
}
/**
* Evaluates the continued fraction.
* <pre>
* a1
* b0 + ------------------
* b1 + a2
* -------------
* b2 + a3
* --------
* b3 + ...
* </pre>
*
* <p>Setting coefficient a<sub>n</sub> to zero will signal the end of the recursive evaluation.
*
* <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
* Both of the first generated terms a and b are used. This fraction evaluation
* can be used when:
* <ul>
* <li>b<sub>0</sub> is not part of a regular series
* <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component
* <li>b<sub>0</sub> is very small and the result is expected to approach zero
* </ul>
*
* <p><b>Usage Note</b>
*
* <p>This method is not functionally identical to calling
* {@link #value(Supplier, double, int)} with the generator configured to provide term
* "b<sub>0</sub>" in the first coefficient. In some cases the computed result from
* the two variations may be different by more than the provided epsilon. The
* convergence of the continued fraction algorithm relies on computing an update
* multiplier applied to the current value. Convergence is faster if the initial value
* is close to the final value. The {@link #value(Supplier, double, int)} method will
* initialise the current value using b<sub>0</sub> and evaluate the continued
* fraction using updates computed from the generated coefficients. This method
* initialises the algorithm using b1 to evaluate part of the continued fraction and
* computes the result as:
*
* <pre>
* a1
* b0 + ------
* part
* </pre>
*
* <p>This is preferred if b<sub>0</sub> is smaller in magnitude than the continued
* fraction component. In particular the evaluation algorithm sets a bound on the
* minimum initial value as {@code 1e-50}. If b<sub>0</sub> is smaller than this value
* then using this method is the preferred evaluation.
*
* @param b0 Coefficient b<sub>0</sub>.
* @param gen Generator of coefficients.
* @param epsilon Maximum relative error allowed.
* @param maxIterations Maximum number of iterations.
* @return the value of the continued fraction.
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number
* of iterations is reached before the expected convergence is achieved.
* @see #value(Supplier,double,int)
*/
public static double value(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) {
// Use the first b coefficient to seed the evaluation of the fraction.
// Coefficient a is used to compute the final result as the numerator term a1.
// The supplied b0 is added to the result.
final Coefficient c = gen.get();
return b0 + c.getA() / evaluate(c.getB(), gen, epsilon, maxIterations);
}
/**
* Evaluates the continued fraction using the modified Lentz algorithm described in
* Thompson and Barnett (1986) Journal of Computational Physics 64, 490-509.
* <pre>
* a1
* b0 + ------------------
* b1 + a2
* -------------
* b2 + a3
* --------
* b3 + ...
* </pre>
*
* <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
* Both of the first generated terms a and b are used.
*
* <p><b>Implementation Note</b>
*
* <p>This method is private and functionally different from
* {@link #value(double, Supplier, double, int)}. The convergence of the algorithm relies on
* computing an update multiplier applied to the current value, initialised as b0. Accuracy
* of the evaluation can be effected if the magnitude of b0 is very different from later
* terms. In particular if initialised as 0 the algorithm will not function and so must
* set b0 to a small non-zero number. The public methods with the leading b0 term
* provide evaluation of the fraction if the term b0 is zero.
*
* @param b0 Coefficient b<sub>0</sub>.
* @param gen Generator of coefficients.
* @param epsilon Maximum relative error allowed.
* @param maxIterations Maximum number of iterations.
* @return the value of the continued fraction.
* @throws ArithmeticException if the algorithm fails to converge or if the maximal number
* of iterations is reached before the expected convergence is achieved.
*/
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);
}
/**
* Returns the value, or if close to zero returns a small epsilon of the same sign.
*
* <p>This method is used in Thompson & Barnett to monitor both the numerator and denominator
* ratios for approaches to zero.
*
* @param value the value
* @return the value (or small epsilon)
*/
private static double updateIfCloseToZero(double value) {
return Math.abs(value) < SMALL ? Math.copySign(SMALL, value) : value;
}
}