    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    import org.apache.commons.math.exception.util.LocalizedFormats;
    023
    024    /**
    025     * Provides a generic means to evaluate continued fractions.  Subclasses simply
    026     * provided the a and b coefficients to evaluate the continued fraction.
    027     *
    028     * <p>
    029     * References:
    030     * <ul>
    031     * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">
    032     * Continued Fraction</a></li>
    033     * </ul>
    034     * </p>
    035     *
    036     * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $
    037     */
    038    public abstract class ContinuedFraction {
    039
    040        /** Maximum allowed numerical error. */
    041        private static final double DEFAULT_EPSILON = 10e-9;
    042
    043        /**
    044         * Default constructor.
    045         */
    046        protected ContinuedFraction() {
    047            super();
    048        }
    049
    050        /**
    051         * Access the n-th a coefficient of the continued fraction.  Since a can be
    052         * a function of the evaluation point, x, that is passed in as well.
    053         * @param n the coefficient index to retrieve.
    054         * @param x the evaluation point.
    055         * @return the n-th a coefficient.
    056         */
    057        protected abstract double getA(int n, double x);
    058
    059        /**
    060         * Access the n-th b coefficient of the continued fraction.  Since b can be
    061         * a function of the evaluation point, x, that is passed in as well.
    062         * @param n the coefficient index to retrieve.
    063         * @param x the evaluation point.
    064         * @return the n-th b coefficient.
    065         */
    066        protected abstract double getB(int n, double x);
    067
    068        /**
    069         * Evaluates the continued fraction at the value x.
    070         * @param x the evaluation point.
    071         * @return the value of the continued fraction evaluated at x.
    072         * @throws MathException if the algorithm fails to converge.
    073         */
    074        public double evaluate(double x) throws MathException {
    075            return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
    076        }
    077
    078        /**
    079         * Evaluates the continued fraction at the value x.
    080         * @param x the evaluation point.
    081         * @param epsilon maximum error allowed.
    082         * @return the value of the continued fraction evaluated at x.
    083         * @throws MathException if the algorithm fails to converge.
    084         */
    085        public double evaluate(double x, double epsilon) throws MathException {
    086            return evaluate(x, epsilon, Integer.MAX_VALUE);
    087        }
    088
    089        /**
    090         * Evaluates the continued fraction at the value x.
    091         * @param x the evaluation point.
    092         * @param maxIterations maximum number of convergents
    093         * @return the value of the continued fraction evaluated at x.
    094         * @throws MathException if the algorithm fails to converge.
    095         */
    096        public double evaluate(double x, int maxIterations) throws MathException {
    097            return evaluate(x, DEFAULT_EPSILON, maxIterations);
    098        }
    099
    100        /**
    101         * <p>
    102         * Evaluates the continued fraction at the value x.
    103         * </p>
    104         *
    105         * <p>
    106         * The implementation of this method is based on equations 14-17 of:
    107         * <ul>
    108         * <li>
    109         *   Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web
    110         *   Resource. <a target="_blank"
    111         *   href="http://mathworld.wolfram.com/ContinuedFraction.html">
    112         *   http://mathworld.wolfram.com/ContinuedFraction.html</a>
    113         * </li>
    114         * </ul>
    115         * The recurrence relationship defined in those equations can result in
    116         * very large intermediate results which can result in numerical overflow.
    117         * As a means to combat these overflow conditions, the intermediate results
    118         * are scaled whenever they threaten to become numerically unstable.</p>
    119         *
    120         * @param x the evaluation point.
    121         * @param epsilon maximum error allowed.
    122         * @param maxIterations maximum number of convergents
    123         * @return the value of the continued fraction evaluated at x.
    124         * @throws MathException if the algorithm fails to converge.
    125         */
    126        public double evaluate(double x, double epsilon, int maxIterations)
    127            throws MathException
    128        {
    129            double p0 = 1.0;
    130            double p1 = getA(0, x);
    131            double q0 = 0.0;
    132            double q1 = 1.0;
    133            double c = p1 / q1;
    134            int n = 0;
    135            double relativeError = Double.MAX_VALUE;
    136            while (n < maxIterations && relativeError > epsilon) {
    137                ++n;
    138                double a = getA(n, x);
    139                double b = getB(n, x);
    140                double p2 = a * p1 + b * p0;
    141                double q2 = a * q1 + b * q0;
    142                boolean infinite = false;
    143                if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
    144                    /*
    145                     * Need to scale. Try successive powers of the larger of a or b
    146                     * up to 5th power. Throw ConvergenceException if one or both
    147                     * of p2, q2 still overflow.
    148                     */
    149                    double scaleFactor = 1d;
    150                    double lastScaleFactor = 1d;
    151                    final int maxPower = 5;
    152                    final double scale = FastMath.max(a,b);
    153                    if (scale <= 0) {  // Can't scale
    154                        throw new ConvergenceException(
    155                                LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
    156                                 x);
    157                    }
    158                    infinite = true;
    159                    for (int i = 0; i < maxPower; i++) {
    160                        lastScaleFactor = scaleFactor;
    161                        scaleFactor *= scale;
    162                        if (a != 0.0 && a > b) {
    163                            p2 = p1 / lastScaleFactor + (b / scaleFactor * p0);
    164                            q2 = q1 / lastScaleFactor + (b / scaleFactor * q0);
    165                        } else if (b != 0) {
    166                            p2 = (a / scaleFactor * p1) + p0 / lastScaleFactor;
    167                            q2 = (a / scaleFactor * q1) + q0 / lastScaleFactor;
    168                        }
    169                        infinite = Double.isInfinite(p2) || Double.isInfinite(q2);
    170                        if (!infinite) {
    171                            break;
    172                        }
    173                    }
    174                }
    175
    176                if (infinite) {
    177                   // Scaling failed
    178                   throw new ConvergenceException(
    179                     LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
    180                      x);
    181                }
    182
    183                double r = p2 / q2;
    184
    185                if (Double.isNaN(r)) {
    186                    throw new ConvergenceException(
    187                      LocalizedFormats.CONTINUED_FRACTION_NAN_DIVERGENCE,
    188                      x);
    189                }
    190                relativeError = FastMath.abs(r / c - 1.0);
    191
    192                // prepare for next iteration
    193                c = p2 / q2;
    194                p0 = p1;
    195                p1 = p2;
    196                q0 = q1;
    197                q1 = q2;
    198            }
    199
    200            if (n >= maxIterations) {
    201                throw new MaxIterationsExceededException(maxIterations,
    202                    LocalizedFormats.NON_CONVERGENT_CONTINUED_FRACTION,
    203                    x);
    204            }
    205
    206            return c;
    207        }
    208    }