001: /*
002: * Copyright 2003-2004 The Apache Software Foundation.
003: *
004: * Licensed under the Apache License, Version 2.0 (the "License");
005: * you may not use this file except in compliance with the License.
006: * You may obtain a copy of the License at
007: *
008: * http://www.apache.org/licenses/LICENSE-2.0
009: *
010: * Unless required by applicable law or agreed to in writing, software
011: * distributed under the License is distributed on an "AS IS" BASIS,
012: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013: * See the License for the specific language governing permissions and
014: * limitations under the License.
015: */
016: package org.apache.commons.math.util;
017:
018: import java.io.Serializable;
019:
020: import org.apache.commons.math.ConvergenceException;
021: import org.apache.commons.math.MathException;
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: 348888 $ $Date: 2005-11-24 23:21:25 -0700 (Thu, 24 Nov 2005) $
036: */
037: public abstract class ContinuedFraction implements Serializable {
038:
039: /** Serialization UID */
040: private static final long serialVersionUID = 1768555336266158242L;
041:
042: /** Maximum allowed numerical error. */
043: private static final double DEFAULT_EPSILON = 10e-9;
044:
045: /**
046: * Default constructor.
047: */
048: protected ContinuedFraction() {
049: super ();
050: }
051:
052: /**
053: * Access the n-th a coefficient of the continued fraction. Since a can be
054: * a function of the evaluation point, x, that is passed in as well.
055: * @param n the coefficient index to retrieve.
056: * @param x the evaluation point.
057: * @return the n-th a coefficient.
058: */
059: protected abstract double getA(int n, double x);
060:
061: /**
062: * Access the n-th b coefficient of the continued fraction. Since b can be
063: * a function of the evaluation point, x, that is passed in as well.
064: * @param n the coefficient index to retrieve.
065: * @param x the evaluation point.
066: * @return the n-th b coefficient.
067: */
068: protected abstract double getB(int n, double x);
069:
070: /**
071: * Evaluates the continued fraction at the value x.
072: * @param x the evaluation point.
073: * @return the value of the continued fraction evaluated at x.
074: * @throws MathException if the algorithm fails to converge.
075: */
076: public double evaluate(double x) throws MathException {
077: return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
078: }
079:
080: /**
081: * Evaluates the continued fraction at the value x.
082: * @param x the evaluation point.
083: * @param epsilon maximum error allowed.
084: * @return the value of the continued fraction evaluated at x.
085: * @throws MathException if the algorithm fails to converge.
086: */
087: public double evaluate(double x, double epsilon)
088: throws MathException {
089: return evaluate(x, epsilon, Integer.MAX_VALUE);
090: }
091:
092: /**
093: * Evaluates the continued fraction at the value x.
094: * @param x the evaluation point.
095: * @param maxIterations maximum number of convergents
096: * @return the value of the continued fraction evaluated at x.
097: * @throws MathException if the algorithm fails to converge.
098: */
099: public double evaluate(double x, int maxIterations)
100: throws MathException {
101: return evaluate(x, DEFAULT_EPSILON, maxIterations);
102: }
103:
104: /**
105: * <p>
106: * Evaluates the continued fraction at the value x.
107: * </p>
108: *
109: * <p>
110: * The implementation of this method is based on equations 14-17 of:
111: * <ul>
112: * <li>
113: * Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web
114: * Resource. <a target="_blank"
115: * href="http://mathworld.wolfram.com/ContinuedFraction.html">
116: * http://mathworld.wolfram.com/ContinuedFraction.html</a>
117: * </li>
118: * </ul>
119: * The recurrence relationship defined in those equations can result in
120: * very large intermediate results which can result in numerical overflow.
121: * As a means to combat these overflow conditions, the intermediate results
122: * are scaled whenever they threaten to become numerically unstable.
123: *
124: * @param x the evaluation point.
125: * @param epsilon maximum error allowed.
126: * @param maxIterations maximum number of convergents
127: * @return the value of the continued fraction evaluated at x.
128: * @throws MathException if the algorithm fails to converge.
129: */
130: public double evaluate(double x, double epsilon, int maxIterations)
131: throws MathException {
132: double p0 = 1.0;
133: double p1 = getA(0, x);
134: double q0 = 0.0;
135: double q1 = 1.0;
136: double c = p1 / q1;
137: int n = 0;
138: double relativeError = Double.MAX_VALUE;
139: while (n < maxIterations && relativeError > epsilon) {
140: ++n;
141: double a = getA(n, x);
142: double b = getB(n, x);
143: double p2 = a * p1 + b * p0;
144: double q2 = a * q1 + b * q0;
145: if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
146: // need to scale
147: if (a != 0.0) {
148: p2 = p1 + (b / a * p0);
149: q2 = q1 + (b / a * q0);
150: } else if (b != 0) {
151: p2 = (a / b * p1) + p0;
152: q2 = (a / b * q1) + q0;
153: } else {
154: // can not scale an convergent is unbounded.
155: throw new ConvergenceException(
156: "Continued fraction convergents diverged to +/- "
157: + "infinity.");
158: }
159: }
160: double r = p2 / q2;
161: relativeError = Math.abs(r / c - 1.0);
162:
163: // prepare for next iteration
164: c = p2 / q2;
165: p0 = p1;
166: p1 = p2;
167: q0 = q1;
168: q1 = q2;
169: }
170:
171: if (n >= maxIterations) {
172: throw new ConvergenceException(
173: "Continued fraction convergents failed to converge.");
174: }
175:
176: return c;
177: }
178: }
|