View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.numbers.fraction;
18  
19  import java.util.function.Supplier;
20  
21  /**
22   * Provides a means to evaluate
23   * <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">generalized continued fractions</a>.
24   *
25   * <p>The continued fraction uses the following form for the numerator ({@code a}) and
26   * denominator ({@code b}) coefficients:
27   * <pre>
28   *              a1
29   * b0 + ------------------
30   *      b1 +      a2
31   *           -------------
32   *           b2 +    a3
33   *                --------
34   *                b3 + ...
35   * </pre>
36   *
37   * <p>A generator of the coefficients must be provided to evaluate the continued fraction.
38   *
39   * <p>The implementation of the fraction evaluation is based on the modified Lentz algorithm
40   * as described on page 508 in:
41   *
42   * <ul>
43   *   <li>
44   *   I. J. Thompson,  A. R. Barnett (1986).
45   *   "Coulomb and Bessel Functions of Complex Arguments and Order."
46   *   Journal of Computational Physics 64, 490-509.
47   *   <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
48   *   https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
49   *   </li>
50   * </ul>
51   *
52   * @see <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">Wikipedia: Generalized continued fraction</a>
53   * @see <a href="https://en.wikipedia.org/wiki/Generalized_continued_fraction">MathWorld: Generalized continued fraction</a>
54   * @since 1.1
55   */
56  public final class GeneralizedContinuedFraction {
57      /**
58       * The value for any number close to zero.
59       *
60       * <p>"The parameter small should be some non-zero number less than typical values of
61       * eps * |b_n|, e.g., 1e-50".
62       */
63      static final double SMALL = 1e-50;
64      /** Default maximum number of iterations. */
65      static final int DEFAULT_ITERATIONS = Integer.MAX_VALUE;
66      /**
67       * Minimum relative error epsilon. Equal to 1 - Math.nextDown(1.0), or 2^-53.
68       *
69       * <p>The epsilon is used to compare the change in the magnitude of the fraction
70       * convergent to 1.0. In theory eps can be 2^-53 reflecting the smallest reduction in
71       * magnitude possible i.e. {@code next = previous * Math.nextDown(1.0)}, or zero
72       * reflecting exact convergence.
73       *
74       * <p>If set to zero then the algorithm requires exact convergence which may not be possible
75       * due to floating point error in the algorithm. For example the golden ratio will not
76       * converge.
77       *
78       * <p>The minimum value will stop the recursive evaluation at the smallest possible
79       * increase or decrease in the convergent.
80       */
81      private static final double MIN_EPSILON = 0x1.0p-53;
82      /** Maximum relative error epsilon. This is configured to prevent incorrect usage. Values
83       * higher than 1.0 invalidate the relative error lower bound of {@code (1 - eps) / 1}.
84       * Set to 0.5 which is a very weak relative error tolerance. */
85      private static final double MAX_EPSILON = 0.5;
86      /** Default low threshold for change in magnitude. Precomputed using MIN_EPSILON.
87       * Equal to 1 - 2^-53. */
88      private static final double DEFAULT_LOW = 1 - MIN_EPSILON;
89      /** Default absolute difference threshold for change in magnitude. Precomputed using MIN_EPSILON.
90       * Equal to {@code 1 / (1 - 2^-53) = 2^-52}. */
91      private static final double DEFAULT_EPS = 0x1.0p-52;
92  
93      /**
94       * Defines the <a href="https://mathworld.wolfram.com/GeneralizedContinuedFraction.html">
95       * {@code n}-th "a" and "b" coefficients</a> of the continued fraction.
96       *
97       * @since 1.1
98       */
99      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 }