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.math3.util; 018 019 import org.apache.commons.math3.exception.ConvergenceException; 020 import org.apache.commons.math3.exception.MaxCountExceededException; 021 import org.apache.commons.math3.exception.util.LocalizedFormats; 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 $Id: ContinuedFraction.java 1416643 2012-12-03 19:37:14Z tn $ 036 */ 037 public abstract class ContinuedFraction { 038 /** Maximum allowed numerical error. */ 039 private static final double DEFAULT_EPSILON = 10e-9; 040 041 /** 042 * Default constructor. 043 */ 044 protected ContinuedFraction() { 045 super(); 046 } 047 048 /** 049 * Access the n-th a coefficient of the continued fraction. Since a can be 050 * a function of the evaluation point, x, that is passed in as well. 051 * @param n the coefficient index to retrieve. 052 * @param x the evaluation point. 053 * @return the n-th a coefficient. 054 */ 055 protected abstract double getA(int n, double x); 056 057 /** 058 * Access the n-th b coefficient of the continued fraction. Since b can be 059 * a function of the evaluation point, x, that is passed in as well. 060 * @param n the coefficient index to retrieve. 061 * @param x the evaluation point. 062 * @return the n-th b coefficient. 063 */ 064 protected abstract double getB(int n, double x); 065 066 /** 067 * Evaluates the continued fraction at the value x. 068 * @param x the evaluation point. 069 * @return the value of the continued fraction evaluated at x. 070 * @throws ConvergenceException if the algorithm fails to converge. 071 */ 072 public double evaluate(double x) throws ConvergenceException { 073 return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE); 074 } 075 076 /** 077 * Evaluates the continued fraction at the value x. 078 * @param x the evaluation point. 079 * @param epsilon maximum error allowed. 080 * @return the value of the continued fraction evaluated at x. 081 * @throws ConvergenceException if the algorithm fails to converge. 082 */ 083 public double evaluate(double x, double epsilon) throws ConvergenceException { 084 return evaluate(x, epsilon, Integer.MAX_VALUE); 085 } 086 087 /** 088 * Evaluates the continued fraction at the value x. 089 * @param x the evaluation point. 090 * @param maxIterations maximum number of convergents 091 * @return the value of the continued fraction evaluated at x. 092 * @throws ConvergenceException if the algorithm fails to converge. 093 * @throws MaxCountExceededException if maximal number of iterations is reached 094 */ 095 public double evaluate(double x, int maxIterations) 096 throws ConvergenceException, MaxCountExceededException { 097 return evaluate(x, DEFAULT_EPSILON, maxIterations); 098 } 099 100 /** 101 * Evaluates the continued fraction at the value x. 102 * <p> 103 * The implementation of this method is based on the modified Lentz algorithm as described 104 * on page 18 ff. in: 105 * <ul> 106 * <li> 107 * I. J. Thompson, A. R. Barnett. "Coulomb and Bessel Functions of Complex Arguments and Order." 108 * <a target="_blank" href="http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf"> 109 * http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a> 110 * </li> 111 * </ul> 112 * <b>Note:</b> the implementation uses the terms a<sub>i</sub> and b<sub>i</sub> as defined in 113 * <a href="http://mathworld.wolfram.com/ContinuedFraction.html">Continued Fraction @ MathWorld</a>. 114 * </p> 115 * 116 * @param x the evaluation point. 117 * @param epsilon maximum error allowed. 118 * @param maxIterations maximum number of convergents 119 * @return the value of the continued fraction evaluated at x. 120 * @throws ConvergenceException if the algorithm fails to converge. 121 * @throws MaxCountExceededException if maximal number of iterations is reached 122 */ 123 public double evaluate(double x, double epsilon, int maxIterations) 124 throws ConvergenceException, MaxCountExceededException { 125 final double small = 1e-50; 126 double hPrev = getA(0, x); 127 128 // use the value of small as epsilon criteria for zero checks 129 if (Precision.equals(hPrev, 0.0, small)) { 130 hPrev = small; 131 } 132 133 int n = 1; 134 double dPrev = 0.0; 135 double cPrev = hPrev; 136 double hN = hPrev; 137 138 while (n < maxIterations) { 139 final double a = getA(n, x); 140 final double b = getB(n, x); 141 142 double dN = a + b * dPrev; 143 if (Precision.equals(dN, 0.0, small)) { 144 dN = small; 145 } 146 double cN = a + b / cPrev; 147 if (Precision.equals(cN, 0.0, small)) { 148 cN = small; 149 } 150 151 dN = 1 / dN; 152 final double deltaN = cN * dN; 153 hN = hPrev * deltaN; 154 155 if (Double.isInfinite(hN)) { 156 throw new ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, 157 x); 158 } 159 if (Double.isNaN(hN)) { 160 throw new ConvergenceException(LocalizedFormats.CONTINUED_FRACTION_NAN_DIVERGENCE, 161 x); 162 } 163 164 if (FastMath.abs(deltaN - 1.0) < epsilon) { 165 break; 166 } 167 168 dPrev = dN; 169 cPrev = cN; 170 hPrev = hN; 171 n++; 172 } 173 174 if (n >= maxIterations) { 175 throw new MaxCountExceededException(LocalizedFormats.NON_CONVERGENT_CONTINUED_FRACTION, 176 maxIterations, x); 177 } 178 179 return hN; 180 } 181 182 }