diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/util/ContinuedFraction.java')
-rw-r--r-- | src/main/java/org/apache/commons/math3/util/ContinuedFraction.java | 181 |
1 files changed, 181 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java b/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java new file mode 100644 index 0000000..69b4964 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/util/ContinuedFraction.java @@ -0,0 +1,181 @@ +/* + * 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.math3.util; + +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.util.LocalizedFormats; + +/** + * Provides a generic means to evaluate continued fractions. Subclasses simply provided the a and b + * coefficients to evaluate the continued fraction. + * + * <p>References: + * + * <ul> + * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">Continued Fraction</a> + * </ul> + */ +public abstract class ContinuedFraction { + /** Maximum allowed numerical error. */ + private static final double DEFAULT_EPSILON = 10e-9; + + /** Default constructor. */ + protected ContinuedFraction() { + super(); + } + + /** + * Access the n-th a coefficient of the continued fraction. Since a can be a function of the + * evaluation point, x, that is passed in as well. + * + * @param n the coefficient index to retrieve. + * @param x the evaluation point. + * @return the n-th a coefficient. + */ + protected abstract double getA(int n, double x); + + /** + * Access the n-th b coefficient of the continued fraction. Since b can be a function of the + * evaluation point, x, that is passed in as well. + * + * @param n the coefficient index to retrieve. + * @param x the evaluation point. + * @return the n-th b coefficient. + */ + protected abstract double getB(int n, double x); + + /** + * Evaluates the continued fraction at the value x. + * + * @param x the evaluation point. + * @return the value of the continued fraction evaluated at x. + * @throws ConvergenceException if the algorithm fails to converge. + */ + public double evaluate(double x) throws ConvergenceException { + return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE); + } + + /** + * Evaluates the continued fraction at the value x. + * + * @param x the evaluation point. + * @param epsilon maximum error allowed. + * @return the value of the continued fraction evaluated at x. + * @throws ConvergenceException if the algorithm fails to converge. + */ + public double evaluate(double x, double epsilon) throws ConvergenceException { + return evaluate(x, epsilon, Integer.MAX_VALUE); + } + + /** + * Evaluates the continued fraction at the value x. + * + * @param x the evaluation point. + * @param maxIterations maximum number of convergents + * @return the value of the continued fraction evaluated at x. + * @throws ConvergenceException if the algorithm fails to converge. + * @throws MaxCountExceededException if maximal number of iterations is reached + */ + public double evaluate(double x, int maxIterations) + throws ConvergenceException, MaxCountExceededException { + return evaluate(x, DEFAULT_EPSILON, maxIterations); + } + + /** + * Evaluates the continued fraction at the value x. + * + * <p>The implementation of this method is based on the modified Lentz algorithm as described on + * page 18 ff. in: + * + * <ul> + * <li>I. J. Thompson, A. R. Barnett. "Coulomb and Bessel Functions of Complex Arguments and + * Order." <a target="_blank" + * href="http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf"> + * http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a> + * </ul> + * + * <b>Note:</b> the implementation uses the terms a<sub>i</sub> and b<sub>i</sub> as defined in + * <a href="http://mathworld.wolfram.com/ContinuedFraction.html">Continued Fraction @ + * MathWorld</a>. + * + * @param x the evaluation point. + * @param epsilon maximum error allowed. + * @param maxIterations maximum number of convergents + * @return the value of the continued fraction evaluated at x. + * @throws ConvergenceException if the algorithm fails to converge. + * @throws MaxCountExceededException if maximal number of iterations is reached + */ + public double evaluate(double x, double epsilon, int maxIterations) + throws ConvergenceException, MaxCountExceededException { + final double small = 1e-50; + double hPrev = getA(0, x); + + // use the value of small as epsilon criteria for zero checks + if (Precision.equals(hPrev, 0.0, small)) { + hPrev = small; + } + + int n = 1; + double dPrev = 0.0; + double cPrev = hPrev; + double hN = hPrev; + + while (n < maxIterations) { + final double a = getA(n, x); + final double b = getB(n, x); + + double dN = a + b * dPrev; + if (Precision.equals(dN, 0.0, small)) { + dN = small; + } + double cN = a + b / cPrev; + if (Precision.equals(cN, 0.0, small)) { + cN = small; + } + + dN = 1 / dN; + final double deltaN = cN * dN; + hN = hPrev * deltaN; + + if (Double.isInfinite(hN)) { + throw new ConvergenceException( + LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, x); + } + if (Double.isNaN(hN)) { + throw new ConvergenceException( + LocalizedFormats.CONTINUED_FRACTION_NAN_DIVERGENCE, x); + } + + if (FastMath.abs(deltaN - 1.0) < epsilon) { + break; + } + + dPrev = dN; + cPrev = cN; + hPrev = hN; + n++; + } + + if (n >= maxIterations) { + throw new MaxCountExceededException( + LocalizedFormats.NON_CONVERGENT_CONTINUED_FRACTION, maxIterations, x); + } + + return hN; + } +} |