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 }