1 /* 2 * Copyright 2003-2004 The Apache Software Foundation. 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 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 // need to scale 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 // can not scale an convergent is unbounded. 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 // prepare for next iteration 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 }