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 }