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    }