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: 348888 $ $Date: 2005-11-24 23:21:25 -0700 (Thu, 24 Nov 2005) $
36 */
37 public abstract class ContinuedFraction implements Serializable {
38
39 /** Serialization UID */
40 private 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 * <p>
104 * Evaluates the continued fraction at the value x.
105 * </p>
106 *
107 * <p>
108 * The implementation of this method is based on equations 14-17 of:
109 * <ul>
110 * <li>
111 * Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web
112 * Resource. <a target="_blank"
113 * href="http://mathworld.wolfram.com/ContinuedFraction.html">
114 * http://mathworld.wolfram.com/ContinuedFraction.html</a>
115 * </li>
116 * </ul>
117 * The recurrence relationship defined in those equations can result in
118 * very large intermediate results which can result in numerical overflow.
119 * As a means to combat these overflow conditions, the intermediate results
120 * are scaled whenever they threaten to become numerically unstable.
121 *
122 * @param x the evaluation point.
123 * @param epsilon maximum error allowed.
124 * @param maxIterations maximum number of convergents
125 * @return the value of the continued fraction evaluated at x.
126 * @throws MathException if the algorithm fails to converge.
127 */
128 public double evaluate(double x, double epsilon, int maxIterations)
129 throws MathException
130 {
131 double p0 = 1.0;
132 double p1 = getA(0, x);
133 double q0 = 0.0;
134 double q1 = 1.0;
135 double c = p1 / q1;
136 int n = 0;
137 double relativeError = Double.MAX_VALUE;
138 while (n < maxIterations && relativeError > epsilon) {
139 ++n;
140 double a = getA(n, x);
141 double b = getB(n, x);
142 double p2 = a * p1 + b * p0;
143 double q2 = a * q1 + b * q0;
144 if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
145
146 if (a != 0.0) {
147 p2 = p1 + (b / a * p0);
148 q2 = q1 + (b / a * q0);
149 } else if (b != 0) {
150 p2 = (a / b * p1) + p0;
151 q2 = (a / b * q1) + q0;
152 } else {
153
154 throw new ConvergenceException(
155 "Continued fraction convergents diverged to +/- " +
156 "infinity.");
157 }
158 }
159 double r = p2 / q2;
160 relativeError = Math.abs(r / c - 1.0);
161
162
163 c = p2 / q2;
164 p0 = p1;
165 p1 = p2;
166 q0 = q1;
167 q1 = q2;
168 }
169
170 if (n >= maxIterations) {
171 throw new ConvergenceException(
172 "Continued fraction convergents failed to converge.");
173 }
174
175 return c;
176 }
177 }