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 import org.apache.commons.math.exception.util.LocalizedFormats; 23 24 /** 25 * Provides a generic means to evaluate continued fractions. Subclasses simply 26 * provided the a and b coefficients to evaluate the continued fraction. 27 * 28 * <p> 29 * References: 30 * <ul> 31 * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html"> 32 * Continued Fraction</a></li> 33 * </ul> 34 * </p> 35 * 36 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $ 37 */ 38 public abstract class ContinuedFraction { 39 40 /** Maximum allowed numerical error. */ 41 private static final double DEFAULT_EPSILON = 10e-9; 42 43 /** 44 * Default constructor. 45 */ 46 protected ContinuedFraction() { 47 super(); 48 } 49 50 /** 51 * Access the n-th a coefficient of the continued fraction. Since a can be 52 * a function of the evaluation point, x, that is passed in as well. 53 * @param n the coefficient index to retrieve. 54 * @param x the evaluation point. 55 * @return the n-th a coefficient. 56 */ 57 protected abstract double getA(int n, double x); 58 59 /** 60 * Access the n-th b coefficient of the continued fraction. Since b can be 61 * a function of the evaluation point, x, that is passed in as well. 62 * @param n the coefficient index to retrieve. 63 * @param x the evaluation point. 64 * @return the n-th b coefficient. 65 */ 66 protected abstract double getB(int n, double x); 67 68 /** 69 * Evaluates the continued fraction at the value x. 70 * @param x the evaluation point. 71 * @return the value of the continued fraction evaluated at x. 72 * @throws MathException if the algorithm fails to converge. 73 */ 74 public double evaluate(double x) throws MathException { 75 return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE); 76 } 77 78 /** 79 * Evaluates the continued fraction at the value x. 80 * @param x the evaluation point. 81 * @param epsilon maximum error allowed. 82 * @return the value of the continued fraction evaluated at x. 83 * @throws MathException if the algorithm fails to converge. 84 */ 85 public double evaluate(double x, double epsilon) throws MathException { 86 return evaluate(x, epsilon, Integer.MAX_VALUE); 87 } 88 89 /** 90 * Evaluates the continued fraction at the value x. 91 * @param x the evaluation point. 92 * @param maxIterations maximum number of convergents 93 * @return the value of the continued fraction evaluated at x. 94 * @throws MathException if the algorithm fails to converge. 95 */ 96 public double evaluate(double x, int maxIterations) throws MathException { 97 return evaluate(x, DEFAULT_EPSILON, maxIterations); 98 } 99 100 /** 101 * <p> 102 * Evaluates the continued fraction at the value x. 103 * </p> 104 * 105 * <p> 106 * The implementation of this method is based on equations 14-17 of: 107 * <ul> 108 * <li> 109 * Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web 110 * Resource. <a target="_blank" 111 * href="http://mathworld.wolfram.com/ContinuedFraction.html"> 112 * http://mathworld.wolfram.com/ContinuedFraction.html</a> 113 * </li> 114 * </ul> 115 * The recurrence relationship defined in those equations can result in 116 * very large intermediate results which can result in numerical overflow. 117 * As a means to combat these overflow conditions, the intermediate results 118 * are scaled whenever they threaten to become numerically unstable.</p> 119 * 120 * @param x the evaluation point. 121 * @param epsilon maximum error allowed. 122 * @param maxIterations maximum number of convergents 123 * @return the value of the continued fraction evaluated at x. 124 * @throws MathException if the algorithm fails to converge. 125 */ 126 public double evaluate(double x, double epsilon, int maxIterations) 127 throws MathException 128 { 129 double p0 = 1.0; 130 double p1 = getA(0, x); 131 double q0 = 0.0; 132 double q1 = 1.0; 133 double c = p1 / q1; 134 int n = 0; 135 double relativeError = Double.MAX_VALUE; 136 while (n < maxIterations && relativeError > epsilon) { 137 ++n; 138 double a = getA(n, x); 139 double b = getB(n, x); 140 double p2 = a * p1 + b * p0; 141 double q2 = a * q1 + b * q0; 142 boolean infinite = false; 143 if (Double.isInfinite(p2) || Double.isInfinite(q2)) { 144 /* 145 * Need to scale. Try successive powers of the larger of a or b 146 * up to 5th power. Throw ConvergenceException if one or both 147 * of p2, q2 still overflow. 148 */ 149 double scaleFactor = 1d; 150 double lastScaleFactor = 1d; 151 final int maxPower = 5; 152 final double scale = FastMath.max(a,b); 153 if (scale <= 0) { // Can't scale 154 throw new ConvergenceException( 155 LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, 156 x); 157 } 158 infinite = true; 159 for (int i = 0; i < maxPower; i++) { 160 lastScaleFactor = scaleFactor; 161 scaleFactor *= scale; 162 if (a != 0.0 && a > b) { 163 p2 = p1 / lastScaleFactor + (b / scaleFactor * p0); 164 q2 = q1 / lastScaleFactor + (b / scaleFactor * q0); 165 } else if (b != 0) { 166 p2 = (a / scaleFactor * p1) + p0 / lastScaleFactor; 167 q2 = (a / scaleFactor * q1) + q0 / lastScaleFactor; 168 } 169 infinite = Double.isInfinite(p2) || Double.isInfinite(q2); 170 if (!infinite) { 171 break; 172 } 173 } 174 } 175 176 if (infinite) { 177 // Scaling failed 178 throw new ConvergenceException( 179 LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, 180 x); 181 } 182 183 double r = p2 / q2; 184 185 if (Double.isNaN(r)) { 186 throw new ConvergenceException( 187 LocalizedFormats.CONTINUED_FRACTION_NAN_DIVERGENCE, 188 x); 189 } 190 relativeError = FastMath.abs(r / c - 1.0); 191 192 // prepare for next iteration 193 c = p2 / q2; 194 p0 = p1; 195 p1 = p2; 196 q0 = q1; 197 q1 = q2; 198 } 199 200 if (n >= maxIterations) { 201 throw new MaxIterationsExceededException(maxIterations, 202 LocalizedFormats.NON_CONVERGENT_CONTINUED_FRACTION, 203 x); 204 } 205 206 return c; 207 } 208 } 209
