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 */
017package org.apache.commons.math3.util;
018
019import org.apache.commons.math3.exception.ConvergenceException;
020import org.apache.commons.math3.exception.MaxCountExceededException;
021import 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 */
037public 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}