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 import org.apache.commons.numbers.fraction.GeneralizedContinuedFraction.Coefficient;
21
22 /**
23 * Provides a generic means to evaluate
24 * <a href="https://mathworld.wolfram.com/ContinuedFraction.html">continued fractions</a>.
25 *
26 * <p>The continued fraction uses the following form for the numerator ({@code a}) and
27 * denominator ({@code b}) coefficients:
28 * <pre>
29 * a1
30 * b0 + ------------------
31 * b1 + a2
32 * -------------
33 * b2 + a3
34 * --------
35 * b3 + ...
36 * </pre>
37 *
38 * <p>Subclasses must provide the {@link #getA(int,double) a} and {@link #getB(int,double) b}
39 * coefficients to evaluate the continued fraction.
40 *
41 * <p>This class allows evaluation of the fraction for a specified evaluation point {@code x};
42 * the point can be used to express the values of the coefficients.
43 * Evaluation of a continued fraction from a generator of the coefficients can be performed using
44 * {@link GeneralizedContinuedFraction}. This may be preferred if the coefficients can be computed
45 * with updates to the previous coefficients.
46 */
47 public abstract class ContinuedFraction {
48 /**
49 * Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
50 * {@code n}-th "a" coefficient</a> of the continued fraction.
51 *
52 * @param n Index of the coefficient to retrieve.
53 * @param x Evaluation point.
54 * @return the coefficient <code>a<sub>n</sub></code>.
55 */
56 protected abstract double getA(int n, double x);
57
58 /**
59 * Defines the <a href="https://mathworld.wolfram.com/ContinuedFraction.html">
60 * {@code n}-th "b" coefficient</a> of the continued fraction.
61 *
62 * @param n Index of the coefficient to retrieve.
63 * @param x Evaluation point.
64 * @return the coefficient <code>b<sub>n</sub></code>.
65 */
66 protected abstract double getB(int n, double x);
67
68 /**
69 * Evaluates the continued fraction.
70 *
71 * @param x the evaluation point.
72 * @param epsilon Maximum relative error allowed.
73 * @return the value of the continued fraction evaluated at {@code x}.
74 * @throws ArithmeticException if the algorithm fails to converge.
75 * @throws ArithmeticException if the maximal number of iterations is reached
76 * before the expected convergence is achieved.
77 *
78 * @see #evaluate(double,double,int)
79 */
80 public double evaluate(double x, double epsilon) {
81 return evaluate(x, epsilon, GeneralizedContinuedFraction.DEFAULT_ITERATIONS);
82 }
83
84 /**
85 * Evaluates the continued fraction.
86 * <p>
87 * The implementation of this method is based on the modified Lentz algorithm as described
88 * on page 508 in:
89 * </p>
90 *
91 * <ul>
92 * <li>
93 * I. J. Thompson, A. R. Barnett (1986).
94 * "Coulomb and Bessel Functions of Complex Arguments and Order."
95 * Journal of Computational Physics 64, 490-509.
96 * <a target="_blank" href="https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
97 * https://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
98 * </li>
99 * </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 }