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.math.util;
18
19 import org.apache.commons.math.ConvergenceException;
20 import org.apache.commons.math.MathException;
21 import org.apache.commons.math.MaxIterationsExceededException;
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: 778522 $ $Date: 2009-05-25 18:24:50 -0400 (Mon, 25 May 2009) $
36 */
37 public abstract class ContinuedFraction {
38
39 /** Maximum allowed numerical error. */
40 private static final double DEFAULT_EPSILON = 10e-9;
41
42 /**
43 * Default constructor.
44 */
45 protected ContinuedFraction() {
46 super();
47 }
48
49 /**
50 * Access the n-th a coefficient of the continued fraction. Since a can be
51 * a function of the evaluation point, x, that is passed in as well.
52 * @param n the coefficient index to retrieve.
53 * @param x the evaluation point.
54 * @return the n-th a coefficient.
55 */
56 protected abstract double getA(int n, double x);
57
58 /**
59 * Access the n-th b coefficient of the continued fraction. Since b can be
60 * a function of the evaluation point, x, that is passed in as well.
61 * @param n the coefficient index to retrieve.
62 * @param x the evaluation point.
63 * @return the n-th b coefficient.
64 */
65 protected abstract double getB(int n, double x);
66
67 /**
68 * Evaluates the continued fraction at the value x.
69 * @param x the evaluation point.
70 * @return the value of the continued fraction evaluated at x.
71 * @throws MathException if the algorithm fails to converge.
72 */
73 public double evaluate(double x) throws MathException {
74 return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
75 }
76
77 /**
78 * Evaluates the continued fraction at the value x.
79 * @param x the evaluation point.
80 * @param epsilon maximum error allowed.
81 * @return the value of the continued fraction evaluated at x.
82 * @throws MathException if the algorithm fails to converge.
83 */
84 public double evaluate(double x, double epsilon) throws MathException {
85 return evaluate(x, epsilon, Integer.MAX_VALUE);
86 }
87
88 /**
89 * Evaluates the continued fraction at the value x.
90 * @param x the evaluation point.
91 * @param maxIterations maximum number of convergents
92 * @return the value of the continued fraction evaluated at x.
93 * @throws MathException if the algorithm fails to converge.
94 */
95 public double evaluate(double x, int maxIterations) throws MathException {
96 return evaluate(x, DEFAULT_EPSILON, maxIterations);
97 }
98
99 /**
100 * <p>
101 * Evaluates the continued fraction at the value x.
102 * </p>
103 *
104 * <p>
105 * The implementation of this method is based on equations 14-17 of:
106 * <ul>
107 * <li>
108 * Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web
109 * Resource. <a target="_blank"
110 * href="http://mathworld.wolfram.com/ContinuedFraction.html">
111 * http://mathworld.wolfram.com/ContinuedFraction.html</a>
112 * </li>
113 * </ul>
114 * The recurrence relationship defined in those equations can result in
115 * very large intermediate results which can result in numerical overflow.
116 * As a means to combat these overflow conditions, the intermediate results
117 * are scaled whenever they threaten to become numerically unstable.</p>
118 *
119 * @param x the evaluation point.
120 * @param epsilon maximum error allowed.
121 * @param maxIterations maximum number of convergents
122 * @return the value of the continued fraction evaluated at x.
123 * @throws MathException if the algorithm fails to converge.
124 */
125 public double evaluate(double x, double epsilon, int maxIterations)
126 throws MathException
127 {
128 double p0 = 1.0;
129 double p1 = getA(0, x);
130 double q0 = 0.0;
131 double q1 = 1.0;
132 double c = p1 / q1;
133 int n = 0;
134 double relativeError = Double.MAX_VALUE;
135 while (n < maxIterations && relativeError > epsilon) {
136 ++n;
137 double a = getA(n, x);
138 double b = getB(n, x);
139 double p2 = a * p1 + b * p0;
140 double q2 = a * q1 + b * q0;
141 if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
142 // need to scale
143 if (a != 0.0) {
144 p2 = p1 + (b / a * p0);
145 q2 = q1 + (b / a * q0);
146 } else if (b != 0) {
147 p2 = (a / b * p1) + p0;
148 q2 = (a / b * q1) + q0;
149 } else {
150 // can not scale an convergent is unbounded.
151 throw new ConvergenceException(
152 "Continued fraction convergents diverged to +/- infinity for value {0}",
153 x);
154 }
155 }
156 double r = p2 / q2;
157 relativeError = Math.abs(r / c - 1.0);
158
159 // prepare for next iteration
160 c = p2 / q2;
161 p0 = p1;
162 p1 = p2;
163 q0 = q1;
164 q1 = q2;
165 }
166
167 if (n >= maxIterations) {
168 throw new MaxIterationsExceededException(maxIterations,
169 "Continued fraction convergents failed to converge for value {0}",
170 x);
171 }
172
173 return c;
174 }
175 }