1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 package org.apache.commons.math.util;
17
18 import java.io.Serializable;
19
20 import org.apache.commons.math.ConvergenceException;
21 import org.apache.commons.math.MathException;
22
23 /***
24 * Provides a generic means to evaluate continued fractions. Subclasses simply
25 * provided the a and b coefficients to evaluate the continued fraction.
26 *
27 * <p>
28 * References:
29 * <ul>
30 * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">
31 * Continued Fraction</a></li>
32 * </ul>
33 * </p>
34 *
35 * @version $Revision: 1.14 $ $Date: 2004/06/23 16:26:16 $
36 */
37 public abstract class ContinuedFraction implements Serializable {
38
39 /*** Serialization UID */
40 static final long serialVersionUID = 1768555336266158242L;
41
42 /*** Maximum allowed numerical error. */
43 private static final double DEFAULT_EPSILON = 10e-9;
44
45 /***
46 * Default constructor.
47 */
48 protected ContinuedFraction() {
49 super();
50 }
51
52 /***
53 * Access the n-th a coefficient of the continued fraction. Since a can be
54 * a function of the evaluation point, x, that is passed in as well.
55 * @param n the coefficient index to retrieve.
56 * @param x the evaluation point.
57 * @return the n-th a coefficient.
58 */
59 protected abstract double getA(int n, double x);
60
61 /***
62 * Access the n-th b coefficient of the continued fraction. Since b can be
63 * a function of the evaluation point, x, that is passed in as well.
64 * @param n the coefficient index to retrieve.
65 * @param x the evaluation point.
66 * @return the n-th b coefficient.
67 */
68 protected abstract double getB(int n, double x);
69
70 /***
71 * Evaluates the continued fraction at the value x.
72 * @param x the evaluation point.
73 * @return the value of the continued fraction evaluated at x.
74 * @throws MathException if the algorithm fails to converge.
75 */
76 public double evaluate(double x) throws MathException {
77 return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
78 }
79
80 /***
81 * Evaluates the continued fraction at the value x.
82 * @param x the evaluation point.
83 * @param epsilon maximum error allowed.
84 * @return the value of the continued fraction evaluated at x.
85 * @throws MathException if the algorithm fails to converge.
86 */
87 public double evaluate(double x, double epsilon) throws MathException {
88 return evaluate(x, epsilon, Integer.MAX_VALUE);
89 }
90
91 /***
92 * Evaluates the continued fraction at the value x.
93 * @param x the evaluation point.
94 * @param maxIterations maximum number of convergents
95 * @return the value of the continued fraction evaluated at x.
96 * @throws MathException if the algorithm fails to converge.
97 */
98 public double evaluate(double x, int maxIterations) throws MathException {
99 return evaluate(x, DEFAULT_EPSILON, maxIterations);
100 }
101
102 /***
103 * Evaluates the continued fraction at the value x.
104 *
105 * The implementation of this method is based on:
106 * <ul>
107 * <li>O. E-gecio-glu, C . K. Koc, J. Rifa i Coma,
108 * <a href="http://citeseer.ist.psu.edu/egecioglu91fast.html">
109 * On Fast Computation of Continued Fractions</a>, Computers Math. Applic.,
110 * 21(2--3), 1991, 167--169.</li>
111 * </ul>
112 *
113 * @param x the evaluation point.
114 * @param epsilon maximum error allowed.
115 * @param maxIterations maximum number of convergents
116 * @return the value of the continued fraction evaluated at x.
117 * @throws MathException if the algorithm fails to converge.
118 */
119 public double evaluate(double x, double epsilon, int maxIterations)
120 throws MathException
121 {
122 double[][] f = new double[2][2];
123 double[][] a = new double[2][2];
124 double[][] an = new double[2][2];
125
126 a[0][0] = getA(0, x);
127 a[0][1] = 1.0;
128 a[1][0] = 1.0;
129 a[1][1] = 0.0;
130
131 return evaluate(1, x, a, an, f, epsilon, maxIterations);
132 }
133
134 /***
135 * Evaluates the n-th convergent, fn = pn / qn, for this continued fraction
136 * at the value x.
137 * @param n the convergent to compute.
138 * @param x the evaluation point.
139 * @param a (n-1)-th convergent matrix. (Input)
140 * @param an the n-th coefficient matrix. (Output)
141 * @param f the n-th convergent matrix. (Output)
142 * @param epsilon maximum error allowed.
143 * @param maxIterations maximum number of convergents
144 * @return the value of the the n-th convergent for this continued fraction
145 * evaluated at x.
146 * @throws MathException if the algorithm fails to converge.
147 */
148 private double evaluate(
149 int n,
150 double x,
151 double[][] a,
152 double[][] an,
153 double[][] f,
154 double epsilon,
155 int maxIterations) throws MathException
156 {
157 double ret;
158
159
160 an[0][0] = getA(n, x);
161 an[0][1] = 1.0;
162 an[1][0] = getB(n, x);
163 an[1][1] = 0.0;
164
165
166 f[0][0] = (a[0][0] * an[0][0]) + (a[0][1] * an[1][0]);
167 f[0][1] = (a[0][0] * an[0][1]) + (a[0][1] * an[1][1]);
168 f[1][0] = (a[1][0] * an[0][0]) + (a[1][1] * an[1][0]);
169 f[1][1] = (a[1][0] * an[0][1]) + (a[1][1] * an[1][1]);
170
171
172 if (Math.abs((f[0][0] * f[1][1]) - (f[1][0] * f[0][1])) <
173 Math.abs(epsilon * f[1][0] * f[1][1]))
174 {
175 ret = f[0][0] / f[1][0];
176 } else {
177 if (n >= maxIterations) {
178 throw new ConvergenceException(
179 "Continued fraction convergents failed to converge.");
180 }
181
182 ret = evaluate(n + 1, x, f
183 , an
184 , a
185 , epsilon, maxIterations);
186 }
187
188 return ret;
189 }
190 }