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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.util;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalStateException;
26
27 /**
28 * Provides a generic means to evaluate continued fractions. Subclasses simply
29 * provided the a and b coefficients to evaluate the continued fraction.
30 * <p>
31 * References:
32 * <ul>
33 * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">
34 * Continued Fraction</a></li>
35 * </ul>
36 */
37 public abstract class ContinuedFraction {
38 /** Maximum allowed numerical error. */
39 private static final double DEFAULT_EPSILON = 10e-9;
40
41 /**
42 * Default constructor.
43 */
44 protected ContinuedFraction() {
45 super();
46 }
47
48 /**
49 * Access the n-th a coefficient of the continued fraction. Since a can be
50 * a function of the evaluation point, x, that is passed in as well.
51 * @param n the coefficient index to retrieve.
52 * @param x the evaluation point.
53 * @return the n-th a coefficient.
54 */
55 protected abstract double getA(int n, double x);
56
57 /**
58 * Access the n-th b coefficient of the continued fraction. Since b can be
59 * a function of the evaluation point, x, that is passed in as well.
60 * @param n the coefficient index to retrieve.
61 * @param x the evaluation point.
62 * @return the n-th b coefficient.
63 */
64 protected abstract double getB(int n, double x);
65
66 /**
67 * Evaluates the continued fraction at the value x.
68 * @param x the evaluation point.
69 * @return the value of the continued fraction evaluated at x.
70 * @throws MathIllegalStateException if the algorithm fails to converge.
71 */
72 public double evaluate(double x) throws MathIllegalStateException {
73 return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
74 }
75
76 /**
77 * Evaluates the continued fraction at the value x.
78 * @param x the evaluation point.
79 * @param epsilon maximum error allowed.
80 * @return the value of the continued fraction evaluated at x.
81 * @throws MathIllegalStateException if the algorithm fails to converge.
82 */
83 public double evaluate(double x, double epsilon) throws MathIllegalStateException {
84 return evaluate(x, epsilon, Integer.MAX_VALUE);
85 }
86
87 /**
88 * Evaluates the continued fraction at the value x.
89 * @param x the evaluation point.
90 * @param maxIterations maximum number of convergents
91 * @return the value of the continued fraction evaluated at x.
92 * @throws MathIllegalStateException if the algorithm fails to converge.
93 * @throws MathIllegalStateException if maximal number of iterations is reached
94 */
95 public double evaluate(double x, int maxIterations)
96 throws MathIllegalStateException {
97 return evaluate(x, DEFAULT_EPSILON, maxIterations);
98 }
99
100 /**
101 * Evaluates the continued fraction at the value x.
102 * <p>
103 * The implementation of this method is based on the modified Lentz algorithm as described
104 * on page 18 ff. in:
105 * </p>
106 * <ul>
107 * <li>
108 * I. J. Thompson, A. R. Barnett. "Coulomb and Bessel Functions of Complex Arguments and Order."
109 * <a target="_blank" href="http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
110 * http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
111 * </li>
112 * </ul>
113 * <p>
114 * <b>Note:</b> the implementation uses the terms a<sub>i</sub> and b<sub>i</sub> as defined in
115 * <a href="http://mathworld.wolfram.com/ContinuedFraction.html">Continued Fraction @ MathWorld</a>.
116 * </p>
117 *
118 * @param x the evaluation point.
119 * @param epsilon maximum error allowed.
120 * @param maxIterations maximum number of convergents
121 * @return the value of the continued fraction evaluated at x.
122 * @throws MathIllegalStateException if the algorithm fails to converge.
123 * @throws MathIllegalStateException if maximal number of iterations is reached
124 */
125 public double evaluate(double x, double epsilon, int maxIterations)
126 throws MathIllegalStateException {
127 final double small = 1e-50;
128 double hPrev = getA(0, x);
129
130 // use the value of small as epsilon criteria for zero checks
131 if (Precision.equals(hPrev, 0.0, small)) {
132 hPrev = small;
133 }
134
135 int n = 1;
136 double dPrev = 0.0;
137 double cPrev = hPrev;
138 double hN = hPrev;
139
140 while (n < maxIterations) {
141 final double a = getA(n, x);
142 final double b = getB(n, x);
143
144 double dN = a + b * dPrev;
145 if (Precision.equals(dN, 0.0, small)) {
146 dN = small;
147 }
148 double cN = a + b / cPrev;
149 if (Precision.equals(cN, 0.0, small)) {
150 cN = small;
151 }
152
153 dN = 1 / dN;
154 final double deltaN = cN * dN;
155 hN = hPrev * deltaN;
156
157 if (Double.isInfinite(hN)) {
158 throw new MathIllegalStateException(LocalizedCoreFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, x);
159 }
160 if (Double.isNaN(hN)) {
161 throw new MathIllegalStateException(LocalizedCoreFormats.CONTINUED_FRACTION_NAN_DIVERGENCE, x);
162 }
163
164 if (FastMath.abs(deltaN - 1.0) < epsilon) {
165 break;
166 }
167
168 dPrev = dN;
169 cPrev = cN;
170 hPrev = hN;
171 n++;
172 }
173
174 if (n >= maxIterations) {
175 throw new MathIllegalStateException(LocalizedCoreFormats.NON_CONVERGENT_CONTINUED_FRACTION,
176 maxIterations, x);
177 }
178
179 return hN;
180 }
181
182 }