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: 778522 $ $Date: 2009-05-25 18:24:50 -0400 (Mon, 25 May 2009) $
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                if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
142                    // need to scale
143                    if (a != 0.0) {
144                        p2 = p1 + (b / a * p0);
145                        q2 = q1 + (b / a * q0);
146                    } else if (b != 0) {
147                        p2 = (a / b * p1) + p0;
148                        q2 = (a / b * q1) + q0;
149                    } else {
150                        // can not scale an convergent is unbounded.
151                        throw new ConvergenceException(
152                            "Continued fraction convergents diverged to +/- infinity for value {0}",
153                            x);
154                    }
155                }
156                double r = p2 / q2;
157                relativeError = Math.abs(r / c - 1.0);
158                    
159                // prepare for next iteration
160                c = p2 / q2;
161                p0 = p1;
162                p1 = p2;
163                q0 = q1;
164                q1 = q2;
165            }
166    
167            if (n >= maxIterations) {
168                throw new MaxIterationsExceededException(maxIterations,
169                    "Continued fraction convergents failed to converge for value {0}",
170                    x);
171            }
172    
173            return c;
174        }
175    }