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 java.util.function.Supplier;
020import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;
021
022/**
023 * Provides a generic means to evaluate
024 * <a href="https://mathworld.wolfram.com/ContinuedFraction.html">continued fractions</a>.
025 *
026 * <p>The continued fraction uses the following form for the numerator ({@code a}) and
027 * denominator ({@code b}) coefficients:
028 * <pre>
029 *              a1
030 * b0 + ------------------
031 *      b1 +      a2
032 *           -------------
033 *           b2 +    a3
034 *                --------
035 *                b3 + ...
036 * </pre>
037 *
038 * <p>Subclasses must provide the {@link #getA(int,double) a} and {@link #getB(int,double) b}
039 * coefficients to evaluate the continued fraction.
040 *
041 * <p>This class allows evaluation of the fraction for a specified evaluation point {@code x};
042 * the point can be used to express the values of the coefficients.
043 * Evaluation of a continued fraction from a generator of the coefficients can be performed using
044 * {@link GeneralizedContinuedFraction}. This may be preferred if the coefficients can be computed
045 * with updates to the previous coefficients.
046 */
047public abstract class ContinuedFraction {
048    /** Create an instance. */
049    public ContinuedFraction() {}
050
051    /**
052     * Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
053     * {@code n}-th "a" coefficient</a> of the continued fraction.
054     *
055     * @param n Index of the coefficient to retrieve.
056     * @param x Evaluation point.
057     * @return the coefficient <code>a<sub>n</sub></code>.
058     */
059    protected abstract double getA(int n, double x);
060
061    /**
062     * Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
063     * {@code n}-th "b" coefficient</a> of the continued fraction.
064     *
065     * @param n Index of the coefficient to retrieve.
066     * @param x Evaluation point.
067     * @return the coefficient <code>b<sub>n</sub></code>.
068     */
069    protected abstract double getB(int n, double x);
070
071    /**
072     * Evaluates the continued fraction.
073     *
074     * @param x the evaluation point.
075     * @param epsilon Maximum relative error allowed.
076     * @return the value of the continued fraction evaluated at {@code x}.
077     * @throws ArithmeticException if the algorithm fails to converge.
078     * @throws ArithmeticException if the maximal number of iterations is reached
079     * before the expected convergence is achieved.
080     *
081     * @see #evaluate(double,double,int)
082     */
083    public double evaluate(double x, double epsilon) {
084        return evaluate(x, epsilon, GeneralizedContinuedFraction.DEFAULT_ITERATIONS);
085    }
086
087    /**
088     * Evaluates the continued fraction.
089     * <p>
090     * The implementation of this method is based on the modified Lentz algorithm as described
091     * on page 508 in:
092     * </p>
093     *
094     * <ul>
095     *   <li>
096     *   I. J. Thompson,  A. R. Barnett (1986).
097     *   "Coulomb and Bessel Functions of Complex Arguments and Order."
098     *   Journal of Computational Physics 64, 490-509.
099     *   <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
100     *   https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
101     *   </li>
102     * </ul>
103     *
104     * @param x Point at which to evaluate the continued fraction.
105     * @param epsilon Maximum relative error allowed.
106     * @param maxIterations Maximum number of iterations.
107     * @return the value of the continued fraction evaluated at {@code x}.
108     * @throws ArithmeticException if the algorithm fails to converge.
109     * @throws ArithmeticException if the maximal number of iterations is reached
110     * before the expected convergence is achieved.
111     */
112    public double evaluate(double x, double epsilon, int maxIterations) {
113        // Delegate to GeneralizedContinuedFraction
114
115        // Get the first coefficient
116        final double b0 = getB(0, x);
117
118        // Generate coefficients from (a1,b1)
119        final Supplier<Coefficient> gen = new Supplier<Coefficient>() {
120            /** Coefficient index. */
121            private int n;
122            @Override
123            public Coefficient get() {
124                n++;
125                final double a = getA(n, x);
126                final double b = getB(n, x);
127                return Coefficient.of(a, b);
128            }
129        };
130
131        // Invoke appropriate method based on magnitude of first term.
132
133        // If b0 is too small or zero it is set to a non-zero small number to allow
134        // magnitude updates. Avoid this by adding b0 at the end if b0 is small.
135        //
136        // This handles the use case of a negligible initial term. If b1 is also small
137        // then the evaluation starting at b0 or b1 may converge poorly.
138        // One solution is to manually compute the convergent until it is not small
139        // and then evaluate the fraction from the next term:
140        // h1 = b0 + a1 / b1
141        // h2 = b0 + a1 / (b1 + a2 / b2)
142        // ...
143        // hn not 'small', start generator at (n+1):
144        // value = GeneralizedContinuedFraction.value(hn, gen)
145        // This solution is not implemented to avoid recursive complexity.
146
147        if (Math.abs(b0) < GeneralizedContinuedFraction.SMALL) {
148            // Updates from initial convergent b1 and computes:
149            // b0 + a1 / [  b1 + a2 / (b2 + ... ) ]
150            return GeneralizedContinuedFraction.value(b0, gen, epsilon, maxIterations);
151        }
152
153        // Use the package-private evaluate method.
154        // Calling GeneralizedContinuedFraction.value(gen, epsilon, maxIterations)
155        // requires the generator to start from (a0,b0) and repeats computation of b0
156        // and wastes computation of a0.
157
158        // Updates from initial convergent b0:
159        // b0 + a1 / (b1 + ... )
160        return GeneralizedContinuedFraction.evaluate(b0, gen, epsilon, maxIterations);
161    }
162}