ContinuedFraction.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;
import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;
/**
* Provides a generic means to evaluate
* <a href="https://mathworld.wolfram.com/ContinuedFraction.html">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>Subclasses must provide the {@link #getA(int,double) a} and {@link #getB(int,double) b}
* coefficients to evaluate the continued fraction.
*
* <p>This class allows evaluation of the fraction for a specified evaluation point {@code x};
* the point can be used to express the values of the coefficients.
* Evaluation of a continued fraction from a generator of the coefficients can be performed using
* {@link GeneralizedContinuedFraction}. This may be preferred if the coefficients can be computed
* with updates to the previous coefficients.
*/
public abstract class ContinuedFraction {
/** Create an instance. */
public ContinuedFraction() {}
/**
* Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
* {@code n}-th "a" coefficient</a> of the continued fraction.
*
* @param n Index of the coefficient to retrieve.
* @param x Evaluation point.
* @return the coefficient <code>a<sub>n</sub></code>.
*/
protected abstract double getA(int n, double x);
/**
* Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
* {@code n}-th "b" coefficient</a> of the continued fraction.
*
* @param n Index of the coefficient to retrieve.
* @param x Evaluation point.
* @return the coefficient <code>b<sub>n</sub></code>.
*/
protected abstract double getB(int n, double x);
/**
* Evaluates the continued fraction.
*
* @param x the evaluation point.
* @param epsilon Maximum relative error allowed.
* @return the value of the continued fraction evaluated at {@code x}.
* @throws ArithmeticException if the algorithm fails to converge.
* @throws ArithmeticException if the maximal number of iterations is reached
* before the expected convergence is achieved.
*
* @see #evaluate(double,double,int)
*/
public double evaluate(double x, double epsilon) {
return evaluate(x, epsilon, GeneralizedContinuedFraction.DEFAULT_ITERATIONS);
}
/**
* Evaluates the continued fraction.
* <p>
* The implementation of this method is based on the modified Lentz algorithm as described
* on page 508 in:
* </p>
*
* <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>
*
* @param x Point at which to evaluate the continued fraction.
* @param epsilon Maximum relative error allowed.
* @param maxIterations Maximum number of iterations.
* @return the value of the continued fraction evaluated at {@code x}.
* @throws ArithmeticException if the algorithm fails to converge.
* @throws ArithmeticException if the maximal number of iterations is reached
* before the expected convergence is achieved.
*/
public double evaluate(double x, double epsilon, int maxIterations) {
// Delegate to GeneralizedContinuedFraction
// Get the first coefficient
final double b0 = getB(0, x);
// Generate coefficients from (a1,b1)
final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
/** Coefficient index. */
private int n;
@Override
public Coefficient get() {
n++;
final double a = getA(n, x);
final double b = getB(n, x);
return Coefficient.of(a, b);
}
};
// Invoke appropriate method based on magnitude of first term.
// If b0 is too small or zero it is set to a non-zero small number to allow
// magnitude updates. Avoid this by adding b0 at the end if b0 is small.
//
// This handles the use case of a negligible initial term. If b1 is also small
// then the evaluation starting at b0 or b1 may converge poorly.
// One solution is to manually compute the convergent until it is not small
// and then evaluate the fraction from the next term:
// h1 = b0 + a1 / b1
// h2 = b0 + a1 / (b1 + a2 / b2)
// ...
// hn not 'small', start generator at (n+1):
// value = GeneralizedContinuedFraction.value(hn, gen)
// This solution is not implemented to avoid recursive complexity.
if (Math.abs(b0) < GeneralizedContinuedFraction.SMALL) {
// Updates from initial convergent b1 and computes:
// b0 + a1 / [ b1 + a2 / (b2 + ... ) ]
return GeneralizedContinuedFraction.value(b0, gen, epsilon, maxIterations);
}
// Use the package-private evaluate method.
// Calling GeneralizedContinuedFraction.value(gen, epsilon, maxIterations)
// requires the generator to start from (a0,b0) and repeats computation of b0
// and wastes computation of a0.
// Updates from initial convergent b0:
// b0 + a1 / (b1 + ... )
return GeneralizedContinuedFraction.evaluate(b0, gen, epsilon, maxIterations);
}
}