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;
020
021/**
022 * Provides a means to evaluate
023 * <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">generalized 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>A generator of the coefficients must be provided to evaluate the continued fraction.
038 *
039 * <p>The implementation of the fraction evaluation is based on the modified Lentz algorithm
040 * as described on page 508 in:
041 *
042 * <ul>
043 *   <li>
044 *   I. J. Thompson,  A. R. Barnett (1986).
045 *   "Coulomb and Bessel Functions of Complex Arguments and Order."
046 *   Journal of Computational Physics 64, 490-509.
047 *   <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
048 *   https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
049 *   </li>
050 * </ul>
051 *
052 * @see <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">Wikipedia: Generalized continued fraction</a>
053 * @see <a href="https://en.wikipedia.org/wiki/Generalized_continued_fraction">MathWorld: Generalized continued fraction</a>
054 * @since 1.1
055 */
056public final class GeneralizedContinuedFraction {
057    /**
058     * The value for any number close to zero.
059     *
060     * <p>"The parameter small should be some non-zero number less than typical values of
061     * eps * |b_n|, e.g., 1e-50".
062     */
063    static final double SMALL = 1e-50;
064    /** Default maximum number of iterations. */
065    static final int DEFAULT_ITERATIONS = Integer.MAX_VALUE;
066    /**
067     * Minimum relative error epsilon. Equal to 1 - Math.nextDown(1.0), or 2^-53.
068     *
069     * <p>The epsilon is used to compare the change in the magnitude of the fraction
070     * convergent to 1.0. In theory eps can be 2^-53 reflecting the smallest reduction in
071     * magnitude possible i.e. {@code next = previous * Math.nextDown(1.0)}, or zero
072     * reflecting exact convergence.
073     *
074     * <p>If set to zero then the algorithm requires exact convergence which may not be possible
075     * due to floating point error in the algorithm. For example the golden ratio will not
076     * converge.
077     *
078     * <p>The minimum value will stop the recursive evaluation at the smallest possible
079     * increase or decrease in the convergent.
080     */
081    private static final double MIN_EPSILON = 0x1.0p-53;
082    /** Maximum relative error epsilon. This is configured to prevent incorrect usage. Values
083     * higher than 1.0 invalidate the relative error lower bound of {@code (1 - eps) / 1}.
084     * Set to 0.5 which is a very weak relative error tolerance. */
085    private static final double MAX_EPSILON = 0.5;
086    /** Default low threshold for change in magnitude. Precomputed using MIN_EPSILON.
087     * Equal to 1 - 2^-53. */
088    private static final double DEFAULT_LOW = 1 - MIN_EPSILON;
089    /** Default absolute difference threshold for change in magnitude. Precomputed using MIN_EPSILON.
090     * Equal to {@code 1 / (1 - 2^-53) = 2^-52}. */
091    private static final double DEFAULT_EPS = 0x1.0p-52;
092
093    /**
094     * Defines the <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">
095     * {@code n}-th "a" and "b" coefficients</a> of the continued fraction.
096     *
097     * @since 1.1
098     */
099    public static final class Coefficient {
100        /** "a" coefficient. */
101        private final double a;
102        /** "b" coefficient. */
103        private final double b;
104
105        /**
106         * @param a "a" coefficient
107         * @param b "b" coefficient
108         */
109        private Coefficient(double a, double b) {
110            this.a = a;
111            this.b = b;
112        }
113
114        /**
115         * Returns the {@code n}-th "a" coefficient of the continued fraction.
116         *
117         * @return the coefficient <code>a<sub>n</sub></code>.
118         */
119        public double getA() {
120            return a;
121        }
122
123        /**
124         * Returns the {@code n}-th "b" coefficient of the continued fraction.
125         *
126         * @return the coefficient <code>b<sub>n</sub></code>.
127         */
128        public double getB() {
129            return b;
130        }
131
132        /**
133         * Create a new coefficient.
134         *
135         * @param a "a" coefficient
136         * @param b "b" coefficient
137         * @return the coefficient
138         */
139        public static Coefficient of(double a, double b) {
140            return new Coefficient(a, b);
141        }
142    }
143
144    /** No instances. */
145    private GeneralizedContinuedFraction() {}
146
147    /**
148     * Evaluates the continued fraction.
149     *
150     * <p>Note: The first generated partial numerator a<sub>0</sub> is discarded.
151     *
152     * @param gen Generator of coefficients.
153     * @return the value of the continued fraction.
154     * @throws ArithmeticException if the algorithm fails to converge or if the maximal number of
155     * iterations is reached before the expected convergence is achieved.
156     * @see #value(Supplier,double,int)
157     */
158    public static double value(Supplier<Coefficient> gen) {
159        return value(gen, MIN_EPSILON, DEFAULT_ITERATIONS);
160    }
161
162    /**
163     * Evaluates the continued fraction.
164     *
165     * <p>Note: The first generated partial numerator a<sub>0</sub> is discarded.
166     *
167     * @param gen Generator of coefficients.
168     * @param epsilon Maximum relative error allowed.
169     * @return the value of the continued fraction.
170     * @throws ArithmeticException if the algorithm fails to converge or if the maximal number of
171     * iterations is reached before the expected convergence is achieved.
172     * @see #value(Supplier,double,int)
173     */
174    public static double value(Supplier<Coefficient> gen, double epsilon) {
175        return value(gen, epsilon, DEFAULT_ITERATIONS);
176    }
177
178    /**
179     * Evaluates the continued fraction.
180     * <pre>
181     *              a1
182     * b0 + ------------------
183     *      b1 +      a2
184     *           -------------
185     *           b2 +    a3
186     *                --------
187     *                b3 + ...
188     * </pre>
189     *
190     * <p>Setting coefficient a<sub>n</sub> to zero will signal the end of the recursive evaluation.
191     *
192     * <p>Note: The first generated partial numerator a<sub>0</sub> is discarded.
193     *
194     * <p><b>Usage Note</b>
195     *
196     * <p>This method is not functionally identical to calling
197     * {@link #value(double, Supplier, double, int)} with the generator configured to
198     * provide coefficients from n=1 and supplying b<sub>0</sub> separately. In some cases
199     * the computed result from the two variations may be different by more than the
200     * provided epsilon. The other method should be used if b<sub>0</sub> is zero or very
201     * small. See the corresponding javadoc for details.
202     *
203     * @param gen Generator of coefficients.
204     * @param epsilon Maximum relative error allowed.
205     * @param maxIterations Maximum number of iterations.
206     * @return the value of the continued fraction.
207     * @throws ArithmeticException if the algorithm fails to converge or if the maximal number of
208     * iterations is reached before the expected convergence is achieved.
209     * @see #value(double, Supplier, double, int)
210     */
211    public static double value(Supplier<Coefficient> gen, double epsilon, int maxIterations) {
212        // Use the first b coefficient to seed the evaluation of the fraction.
213        // Coefficient a is discarded.
214        final Coefficient c = gen.get();
215        return evaluate(c.getB(), gen, epsilon, maxIterations);
216    }
217
218    /**
219     * Evaluates the continued fraction.
220     *
221     * <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
222     * Both of the first generated terms a and b are used. This fraction evaluation
223     * can be used when:
224     * <ul>
225     *  <li>b<sub>0</sub> is not part of a regular series
226     *  <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component
227     *  <li>b<sub>0</sub> is very small and the result is expected to approach zero
228     * </ul>
229     *
230     * @param b0 Coefficient b<sub>0</sub>.
231     * @param gen Generator of coefficients.
232     * @return the value of the continued fraction.
233     * @throws ArithmeticException if the algorithm fails to converge or if the maximal number
234     * of iterations is reached before the expected convergence is achieved.
235     * @see #value(double,Supplier,double,int)
236     */
237    public static double value(double b0, Supplier<Coefficient> gen) {
238        return value(b0, gen, MIN_EPSILON, DEFAULT_ITERATIONS);
239    }
240
241    /**
242     * Evaluates the continued fraction.
243     *
244     * <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
245     * Both of the first generated terms a and b are used. This fraction evaluation
246     * can be used when:
247     * <ul>
248     *  <li>b<sub>0</sub> is not part of a regular series
249     *  <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component
250     *  <li>b<sub>0</sub> is very small and the result is expected to approach zero
251     * </ul>
252     *
253     * @param b0 Coefficient b<sub>0</sub>.
254     * @param gen Generator of coefficients.
255     * @param epsilon Maximum relative error allowed.
256     * @return the value of the continued fraction.
257     * @throws ArithmeticException if the algorithm fails to converge or if the maximal number
258     * of iterations is reached before the expected convergence is achieved.
259     * @see #value(double,Supplier,double,int)
260     */
261    public static double value(double b0, Supplier<Coefficient> gen, double epsilon) {
262        return value(b0, gen, epsilon, DEFAULT_ITERATIONS);
263    }
264
265    /**
266     * Evaluates the continued fraction.
267     * <pre>
268     *              a1
269     * b0 + ------------------
270     *      b1 +      a2
271     *           -------------
272     *           b2 +    a3
273     *                --------
274     *                b3 + ...
275     * </pre>
276     *
277     * <p>Setting coefficient a<sub>n</sub> to zero will signal the end of the recursive evaluation.
278     *
279     * <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
280     * Both of the first generated terms a and b are used. This fraction evaluation
281     * can be used when:
282     * <ul>
283     *  <li>b<sub>0</sub> is not part of a regular series
284     *  <li>b<sub>0</sub> is zero and the result will evaluate only the continued fraction component
285     *  <li>b<sub>0</sub> is very small and the result is expected to approach zero
286     * </ul>
287     *
288     * <p><b>Usage Note</b>
289     *
290     * <p>This method is not functionally identical to calling
291     * {@link #value(Supplier, double, int)} with the generator configured to provide term
292     * "b<sub>0</sub>" in the first coefficient. In some cases the computed result from
293     * the two variations may be different by more than the provided epsilon. The
294     * convergence of the continued fraction algorithm relies on computing an update
295     * multiplier applied to the current value. Convergence is faster if the initial value
296     * is close to the final value. The {@link #value(Supplier, double, int)} method will
297     * initialise the current value using b<sub>0</sub> and evaluate the continued
298     * fraction using updates computed from the generated coefficients. This method
299     * initialises the algorithm using b1 to evaluate part of the continued fraction and
300     * computes the result as:
301     *
302     * <pre>
303     *        a1
304     * b0 + ------
305     *       part
306     * </pre>
307     *
308     * <p>This is preferred if b<sub>0</sub> is smaller in magnitude than the continued
309     * fraction component. In particular the evaluation algorithm sets a bound on the
310     * minimum initial value as {@code 1e-50}. If b<sub>0</sub> is smaller than this value
311     * then using this method is the preferred evaluation.
312     *
313     * @param b0 Coefficient b<sub>0</sub>.
314     * @param gen Generator of coefficients.
315     * @param epsilon Maximum relative error allowed.
316     * @param maxIterations Maximum number of iterations.
317     * @return the value of the continued fraction.
318     * @throws ArithmeticException if the algorithm fails to converge or if the maximal number
319     * of iterations is reached before the expected convergence is achieved.
320     * @see #value(Supplier,double,int)
321     */
322    public static double value(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) {
323        // Use the first b coefficient to seed the evaluation of the fraction.
324        // Coefficient a is used to compute the final result as the numerator term a1.
325        // The supplied b0 is added to the result.
326        final Coefficient c = gen.get();
327        return b0 + c.getA() / evaluate(c.getB(), gen, epsilon, maxIterations);
328    }
329
330    /**
331     * Evaluates the continued fraction using the modified Lentz algorithm described in
332     * Thompson and Barnett (1986) Journal of Computational Physics 64, 490-509.
333     * <pre>
334     *              a1
335     * b0 + ------------------
336     *      b1 +      a2
337     *           -------------
338     *           b2 +    a3
339     *                --------
340     *                b3 + ...
341     * </pre>
342     *
343     * <p>Note: The initial term b<sub>0</sub> is supplied as an argument.
344     * Both of the first generated terms a and b are used.
345     *
346     * <p><b>Implementation Note</b>
347     *
348     * <p>This method is private and functionally different from
349     * {@link #value(double, Supplier, double, int)}. The convergence of the algorithm relies on
350     * computing an update multiplier applied to the current value, initialised as b0. Accuracy
351     * of the evaluation can be effected if the magnitude of b0 is very different from later
352     * terms. In particular if initialised as 0 the algorithm will not function and so must
353     * set b0 to a small non-zero number. The public methods with the leading b0 term
354     * provide evaluation of the fraction if the term b0 is zero.
355     *
356     * @param b0 Coefficient b<sub>0</sub>.
357     * @param gen Generator of coefficients.
358     * @param epsilon Maximum relative error allowed.
359     * @param maxIterations Maximum number of iterations.
360     * @return the value of the continued fraction.
361     * @throws ArithmeticException if the algorithm fails to converge or if the maximal number
362     * of iterations is reached before the expected convergence is achieved.
363     */
364    static double evaluate(double b0, Supplier<Coefficient> gen, double epsilon, int maxIterations) {
365        // Relative error epsilon should not be zero to prevent drift in the event
366        // that the update ratio never achieves 1.0.
367
368        // Epsilon is the relative change allowed from 1. Configure the absolute limits so
369        // convergence requires: low <= deltaN <= high
370        // low = 1 - eps
371        // high = 1 / (1 - eps)
372        // High is always further from 1 than low in absolute distance. Do not store high
373        // but store the maximum absolute deviation from 1 for convergence = high - 1.
374        // If this is achieved a second check is made against low.
375        double low;
376        double eps;
377        if (epsilon > MIN_EPSILON && epsilon <= MAX_EPSILON) {
378            low = 1 - epsilon;
379            eps = 1 / low - 1;
380        } else {
381            // Precomputed defaults. Used when epsilon <= MIN_EPSILON
382            low = DEFAULT_LOW;
383            eps = DEFAULT_EPS;
384        }
385
386        double hPrev = updateIfCloseToZero(b0);
387
388        // Notes from Thompson and Barnett:
389        //
390        // Fraction convergent: hn = An / Bn
391        // A(-1) = 1, A0 = b0, B(-1) = 0, B0 = 1
392
393        // Compute the ratios:
394        // Dn = B(n-1) / Bn  = 1 / (an * D(n-1) + bn)
395        // Cn = An / A(n-1)  = an / C(n-1) + bn
396        //
397        // Ratio of successive convergents:
398        // delta n = hn / h(n-1)
399        //         = Cn / Dn
400
401        // Avoid divisors being zero (less than machine precision) by shifting them to e.g. 1e-50.
402
403        double dPrev = 0.0;
404        double cPrev = hPrev;
405
406        for (int n = maxIterations; n > 0; n--) {
407            final Coefficient c = gen.get();
408            final double a = c.getA();
409            final double b = c.getB();
410
411            double dN = updateIfCloseToZero(b + a * dPrev);
412            final double cN = updateIfCloseToZero(b + a / cPrev);
413
414            dN = 1 / dN;
415            final double deltaN = cN * dN;
416            final double hN = hPrev * deltaN;
417
418            // If the fraction is convergent then deltaN -> 1.
419            // Computation of deltaN = 0 or deltaN = big will result in zero or overflow.
420            // Directly check for overflow on hN (this ensures the result is finite).
421
422            if (!Double.isFinite(hN)) {
423                throw new FractionException("Continued fraction diverged to " + hN);
424            }
425
426            // Check for underflow on deltaN. This allows fractions to compute zero
427            // if this is the convergent limit.
428            // Note: deltaN is only zero if dN > 1e-50 / min_value, or 2.02e273.
429            // Since dN is the ratio of convergent denominators this magnitude of
430            // ratio is a presumed to be an error.
431            if (deltaN == 0) {
432                throw new FractionException("Ratio of successive convergents is zero");
433            }
434
435            // Update from Thompson and Barnett to use <= eps in place of < eps.
436            // eps = high - 1
437            // A second check is made to ensure:
438            // low <= deltaN <= high
439            if (Math.abs(deltaN - 1) <= eps && deltaN >= low) {
440                return hN;
441            }
442
443            dPrev = dN;
444            cPrev = cN;
445            hPrev = hN;
446        }
447
448        throw new FractionException("Maximum iterations (%d) exceeded", maxIterations);
449    }
450
451    /**
452     * Returns the value, or if close to zero returns a small epsilon of the same sign.
453     *
454     * <p>This method is used in Thompson &amp; Barnett to monitor both the numerator and denominator
455     * ratios for approaches to zero.
456     *
457     * @param value the value
458     * @return the value (or small epsilon)
459     */
460    private static double updateIfCloseToZero(double value) {
461        return Math.abs(value) < SMALL ? Math.copySign(SMALL, value) : value;
462    }
463}