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