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 & 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 }