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="https://mathworld.wolfram.com/ContinuedFraction.html">continued fractions</a>.
024 *
025 * <p>The continued fraction uses the following form for the numerator ({@code a}) and
026 * denominator ({@code b}) coefficients:
027 * <pre>
028 *              a1
029 * b0 + ------------------
030 *      b1 +      a2
031 *           -------------
032 *           b2 +    a3
033 *                --------
034 *                b3 + ...
035 * </pre>
036 *
037 * <p>Subclasses must provide the {@link #getA(int,double) a} and {@link #getB(int,double) b}
038 * coefficients to evaluate the continued fraction.
039 */
040public abstract class ContinuedFraction {
041    /**
042     * The value for any number close to zero.
043     *
044     * <p>"The parameter small should be some non-zero number less than typical values of
045     * eps * |b_n|, e.g., 1e-50".
046     */
047    private static final double SMALL = 1e-50;
048
049    /**
050     * Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
051     * {@code n}-th "a" coefficient</a> of the continued fraction.
052     *
053     * @param n Index of the coefficient to retrieve.
054     * @param x Evaluation point.
055     * @return the coefficient <code>a<sub>n</sub></code>.
056     */
057    protected abstract double getA(int n, double x);
058
059    /**
060     * Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
061     * {@code n}-th "b" coefficient</a> of the continued fraction.
062     *
063     * @param n Index of the coefficient to retrieve.
064     * @param x Evaluation point.
065     * @return the coefficient <code>b<sub>n</sub></code>.
066     */
067    protected abstract double getB(int n, double x);
068
069    /**
070     * Evaluates the continued fraction.
071     *
072     * @param x the evaluation point.
073     * @param epsilon Maximum error allowed.
074     * @return the value of the continued fraction evaluated at {@code x}.
075     * @throws ArithmeticException if the algorithm fails to converge.
076     * @throws ArithmeticException if the maximal number of iterations is reached
077     * before the expected convergence is achieved.
078     *
079     * @see #evaluate(double,double,int)
080     */
081    public double evaluate(double x, double epsilon) {
082        return evaluate(x, epsilon, Integer.MAX_VALUE);
083    }
084
085    /**
086     * Evaluates the continued fraction.
087     * <p>
088     * The implementation of this method is based on the modified Lentz algorithm as described
089     * on page 508 in:
090     * </p>
091     *
092     * <ul>
093     *   <li>
094     *   I. J. Thompson,  A. R. Barnett (1986).
095     *   "Coulomb and Bessel Functions of Complex Arguments and Order."
096     *   Journal of Computational Physics 64, 490-509.
097     *   <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
098     *   https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
099     *   </li>
100     * </ul>
101     *
102     * @param x Point at which to evaluate the continued fraction.
103     * @param epsilon Maximum error allowed.
104     * @param maxIterations Maximum number of iterations.
105     * @return the value of the continued fraction evaluated at {@code x}.
106     * @throws ArithmeticException if the algorithm fails to converge.
107     * @throws ArithmeticException if the maximal number of iterations is reached
108     * before the expected convergence is achieved.
109     */
110    public double evaluate(double x, double epsilon, int maxIterations) {
111        double hPrev = updateIfCloseToZero(getB(0, x));
112
113        int n = 1;
114        double dPrev = 0.0;
115        double cPrev = hPrev;
116        double hN;
117
118        while (n <= maxIterations) {
119            final double a = getA(n, x);
120            final double b = getB(n, x);
121
122            double dN = updateIfCloseToZero(b + a * dPrev);
123            final double cN = updateIfCloseToZero(b + a / cPrev);
124
125            dN = 1 / dN;
126            final double deltaN = cN * dN;
127            hN = hPrev * deltaN;
128
129            if (Double.isInfinite(hN)) {
130                throw new FractionException(
131                    "Continued fraction convergents diverged to +/- infinity for value {0}", x);
132            }
133            if (Double.isNaN(hN)) {
134                throw new FractionException(
135                    "Continued fraction diverged to NaN for value {0}", x);
136            }
137
138            if (Math.abs(deltaN - 1) < epsilon) {
139                return hN;
140            }
141
142            dPrev = dN;
143            cPrev = cN;
144            hPrev = hN;
145            ++n;
146        }
147
148        throw new FractionException("maximal count ({0}) exceeded", maxIterations);
149    }
150
151    /**
152     * Returns the value, or if close to zero returns a small epsilon.
153     *
154     * <p>This method is used in Thompson & Barnett to monitor both the numerator and denominator
155     * ratios for approaches to zero.
156     *
157     * @param value the value
158     * @return the value (or small epsilon)
159     */
160    private static double updateIfCloseToZero(double value) {
161        return Precision.equals(value, 0.0, SMALL) ? SMALL : value;
162    }
163}