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 }