001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.fraction;
018
019import org.apache.commons.numbers.core.Precision;
020
021/**
022 * Provides a generic means to evaluate
023 * <a href="http://mathworld.wolfram.com/ContinuedFraction.html">continued fractions</a>.
024 * Subclasses must provide the {@link #getA(int,double) a} and {@link #getB(int,double) b}
025 * coefficients to evaluate the continued fraction.
026 */
027public abstract class ContinuedFraction {
028    /**
029     * Defines the <a href="http://mathworld.wolfram.com/ContinuedFraction.html">
030     * {@code n}-th "a" coefficient</a> of the continued fraction.
031     *
032     * @param n Index of the coefficient to retrieve.
033     * @param x Evaluation point.
034     * @return the coefficient <code>a<sub>n</sub></code>.
035     */
036    protected abstract double getA(int n, double x);
037
038    /**
039     * Defines the <a href="http://mathworld.wolfram.com/ContinuedFraction.html">
040     * {@code n}-th "b" coefficient</a> of the continued fraction.
041     *
042     * @param n Index of the coefficient to retrieve.
043     * @param x Evaluation point.
044     * @return the coefficient <code>b<sub>n</sub></code>.
045     */
046    protected abstract double getB(int n, double x);
047
048    /**
049     * Evaluates the continued fraction.
050     *
051     * @param x the evaluation point.
052     * @param epsilon Maximum error allowed.
053     * @return the value of the continued fraction evaluated at {@code x}.
054     * @throws ArithmeticException if the algorithm fails to converge.
055     * @throws ArithmeticException if the maximal number of iterations is reached
056     * before the expected convergence is achieved.
057     *
058     * @see #evaluate(double,double,int)
059     */
060    public double evaluate(double x, double epsilon) {
061        return evaluate(x, epsilon, Integer.MAX_VALUE);
062    }
063
064    /**
065     * Evaluates the continued fraction.
066     * <p>
067     * The implementation of this method is based on the modified Lentz algorithm as described
068     * on page 18 ff. in:
069     * </p>
070     *
071     * <ul>
072     *   <li>
073     *   I. J. Thompson,  A. R. Barnett. "Coulomb and Bessel Functions of Complex Arguments and Order."
074     *   <a target="_blank" href="http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
075     *   http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
076     *   </li>
077     * </ul>
078     *
079     * @param x Point at which to evaluate the continued fraction.
080     * @param epsilon Maximum error allowed.
081     * @param maxIterations Maximum number of iterations.
082     * @return the value of the continued fraction evaluated at {@code x}.
083     * @throws ArithmeticException if the algorithm fails to converge.
084     * @throws ArithmeticException if the maximal number of iterations is reached
085     * before the expected convergence is achieved.
086     */
087    public double evaluate(double x, double epsilon, int maxIterations) {
088        final double small = 1e-50;
089        double hPrev = getA(0, x);
090
091        // use the value of small as epsilon criteria for zero checks
092        if (Precision.equals(hPrev, 0.0, small)) {
093            hPrev = small;
094        }
095
096        int n = 1;
097        double dPrev = 0.0;
098        double cPrev = hPrev;
099        double hN;
100
101        while (n <= maxIterations) {
102            final double a = getA(n, x);
103            final double b = getB(n, x);
104
105            double dN = a + b * dPrev;
106            if (Precision.equals(dN, 0.0, small)) {
107                dN = small;
108            }
109            double cN = a + b / cPrev;
110            if (Precision.equals(cN, 0.0, small)) {
111                cN = small;
112            }
113
114            dN = 1 / dN;
115            final double deltaN = cN * dN;
116            hN = hPrev * deltaN;
117
118            if (Double.isInfinite(hN)) {
119                throw new FractionException("Continued fraction convergents diverged to +/- infinity for value {0}",
120                                               x);
121            }
122            if (Double.isNaN(hN)) {
123                throw new FractionException("Continued fraction diverged to NaN for value {0}",
124                                               x);
125            }
126
127            if (Math.abs(deltaN - 1) < epsilon) {
128                return hN;
129            }
130
131            dPrev = dN;
132            cPrev = cN;
133            hPrev = hN;
134            ++n;
135        }
136
137        throw new FractionException("maximal count ({0}) exceeded", maxIterations);
138    }
139}