    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    }