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