001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.util;
018
019 import org.apache.commons.math.ConvergenceException;
020 import org.apache.commons.math.MathException;
021 import org.apache.commons.math.MaxIterationsExceededException;
022
023 /**
024 * Provides a generic means to evaluate continued fractions. Subclasses simply
025 * provided the a and b coefficients to evaluate the continued fraction.
026 *
027 * <p>
028 * References:
029 * <ul>
030 * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">
031 * Continued Fraction</a></li>
032 * </ul>
033 * </p>
034 *
035 * @version $Revision: 920558 $ $Date: 2010-03-08 17:57:32 -0500 (Mon, 08 Mar 2010) $
036 */
037 public abstract class ContinuedFraction {
038
039 /** Maximum allowed numerical error. */
040 private static final double DEFAULT_EPSILON = 10e-9;
041
042 /**
043 * Default constructor.
044 */
045 protected ContinuedFraction() {
046 super();
047 }
048
049 /**
050 * Access the n-th a coefficient of the continued fraction. Since a can be
051 * a function of the evaluation point, x, that is passed in as well.
052 * @param n the coefficient index to retrieve.
053 * @param x the evaluation point.
054 * @return the n-th a coefficient.
055 */
056 protected abstract double getA(int n, double x);
057
058 /**
059 * Access the n-th b coefficient of the continued fraction. Since b can be
060 * a function of the evaluation point, x, that is passed in as well.
061 * @param n the coefficient index to retrieve.
062 * @param x the evaluation point.
063 * @return the n-th b coefficient.
064 */
065 protected abstract double getB(int n, double x);
066
067 /**
068 * Evaluates the continued fraction at the value x.
069 * @param x the evaluation point.
070 * @return the value of the continued fraction evaluated at x.
071 * @throws MathException if the algorithm fails to converge.
072 */
073 public double evaluate(double x) throws MathException {
074 return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
075 }
076
077 /**
078 * Evaluates the continued fraction at the value x.
079 * @param x the evaluation point.
080 * @param epsilon maximum error allowed.
081 * @return the value of the continued fraction evaluated at x.
082 * @throws MathException if the algorithm fails to converge.
083 */
084 public double evaluate(double x, double epsilon) throws MathException {
085 return evaluate(x, epsilon, Integer.MAX_VALUE);
086 }
087
088 /**
089 * Evaluates the continued fraction at the value x.
090 * @param x the evaluation point.
091 * @param maxIterations maximum number of convergents
092 * @return the value of the continued fraction evaluated at x.
093 * @throws MathException if the algorithm fails to converge.
094 */
095 public double evaluate(double x, int maxIterations) throws MathException {
096 return evaluate(x, DEFAULT_EPSILON, maxIterations);
097 }
098
099 /**
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 boolean infinite = false;
142 if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
143 /*
144 * Need to scale. Try successive powers of the larger of a or b
145 * up to 5th power. Throw ConvergenceException if one or both
146 * of p2, q2 still overflow.
147 */
148 double scaleFactor = 1d;
149 double lastScaleFactor = 1d;
150 final int maxPower = 5;
151 final double scale = Math.max(a,b);
152 if (scale <= 0) { // Can't scale
153 throw new ConvergenceException(
154 "Continued fraction convergents diverged to +/- infinity for value {0}",
155 x);
156 }
157 infinite = true;
158 for (int i = 0; i < maxPower; i++) {
159 lastScaleFactor = scaleFactor;
160 scaleFactor *= scale;
161 if (a != 0.0 && a > b) {
162 p2 = p1 / lastScaleFactor + (b / scaleFactor * p0);
163 q2 = q1 / lastScaleFactor + (b / scaleFactor * q0);
164 } else if (b != 0) {
165 p2 = (a / scaleFactor * p1) + p0 / lastScaleFactor;
166 q2 = (a / scaleFactor * q1) + q0 / lastScaleFactor;
167 }
168 infinite = Double.isInfinite(p2) || Double.isInfinite(q2);
169 if (!infinite) {
170 break;
171 }
172 }
173 }
174
175 if (infinite) {
176 // Scaling failed
177 throw new ConvergenceException(
178 "Continued fraction convergents diverged to +/- infinity for value {0}",
179 x);
180 }
181
182 double r = p2 / q2;
183
184 if (Double.isNaN(r)) {
185 throw new ConvergenceException(
186 "Continued fraction diverged to NaN for value {0}",
187 x);
188 }
189 relativeError = Math.abs(r / c - 1.0);
190
191 // prepare for next iteration
192 c = p2 / q2;
193 p0 = p1;
194 p1 = p2;
195 q0 = q1;
196 q1 = q2;
197 }
198
199 if (n >= maxIterations) {
200 throw new MaxIterationsExceededException(maxIterations,
201 "Continued fraction convergents failed to converge for value {0}",
202 x);
203 }
204
205 return c;
206 }
207 }