View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math.util;
18  
19  import java.io.Serializable;
20  
21  import org.apache.commons.math.ConvergenceException;
22  import org.apache.commons.math.MathException;
23  import org.apache.commons.math.MaxIterationsExceededException;
24  
25  /**
26   * Provides a generic means to evaluate continued fractions.  Subclasses simply
27   * provided the a and b coefficients to evaluate the continued fraction.
28   *
29   * <p>
30   * References:
31   * <ul>
32   * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">
33   * Continued Fraction</a></li>
34   * </ul>
35   * </p>
36   *
37   * @version $Revision: 617953 $ $Date: 2008-02-02 22:54:00 -0700 (Sat, 02 Feb 2008) $
38   */
39  public abstract class ContinuedFraction implements Serializable {
40      
41      /** Serialization UID */
42      private static final long serialVersionUID = 1768555336266158242L;
43      
44      /** Maximum allowed numerical error. */
45      private static final double DEFAULT_EPSILON = 10e-9;
46  
47      /**
48       * Default constructor.
49       */
50      protected ContinuedFraction() {
51          super();
52      }
53  
54      /**
55       * Access the n-th a coefficient of the continued fraction.  Since a can be
56       * a function of the evaluation point, x, that is passed in as well.
57       * @param n the coefficient index to retrieve.
58       * @param x the evaluation point.
59       * @return the n-th a coefficient.
60       */
61      protected abstract double getA(int n, double x);
62  
63      /**
64       * Access the n-th b coefficient of the continued fraction.  Since b can be
65       * a function of the evaluation point, x, that is passed in as well.
66       * @param n the coefficient index to retrieve.
67       * @param x the evaluation point.
68       * @return the n-th b coefficient.
69       */
70      protected abstract double getB(int n, double x);
71  
72      /**
73       * Evaluates the continued fraction at the value x.
74       * @param x the evaluation point.
75       * @return the value of the continued fraction evaluated at x. 
76       * @throws MathException if the algorithm fails to converge.
77       */
78      public double evaluate(double x) throws MathException {
79          return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
80      }
81  
82      /**
83       * Evaluates the continued fraction at the value x.
84       * @param x the evaluation point.
85       * @param epsilon maximum error allowed.
86       * @return the value of the continued fraction evaluated at x. 
87       * @throws MathException if the algorithm fails to converge.
88       */
89      public double evaluate(double x, double epsilon) throws MathException {
90          return evaluate(x, epsilon, Integer.MAX_VALUE);
91      }
92  
93      /**
94       * Evaluates the continued fraction at the value x.
95       * @param x the evaluation point.
96       * @param maxIterations maximum number of convergents
97       * @return the value of the continued fraction evaluated at x. 
98       * @throws MathException if the algorithm fails to converge.
99       */
100     public double evaluate(double x, int maxIterations) throws MathException {
101         return evaluate(x, DEFAULT_EPSILON, maxIterations);
102     }
103 
104     /**
105      * <p>
106      * Evaluates the continued fraction at the value x.
107      * </p>
108      * 
109      * <p>
110      * The implementation of this method is based on equations 14-17 of:
111      * <ul>
112      * <li>
113      *   Eric W. Weisstein. "Continued Fraction." From MathWorld--A Wolfram Web
114      *   Resource. <a target="_blank"
115      *   href="http://mathworld.wolfram.com/ContinuedFraction.html">
116      *   http://mathworld.wolfram.com/ContinuedFraction.html</a>
117      * </li>
118      * </ul>
119      * The recurrence relationship defined in those equations can result in
120      * very large intermediate results which can result in numerical overflow.
121      * As a means to combat these overflow conditions, the intermediate results
122      * are scaled whenever they threaten to become numerically unstable.</p>
123      *   
124      * @param x the evaluation point.
125      * @param epsilon maximum error allowed.
126      * @param maxIterations maximum number of convergents
127      * @return the value of the continued fraction evaluated at x. 
128      * @throws MathException if the algorithm fails to converge.
129      */
130     public double evaluate(double x, double epsilon, int maxIterations)
131         throws MathException
132     {
133         double p0 = 1.0;
134         double p1 = getA(0, x);
135         double q0 = 0.0;
136         double q1 = 1.0;
137         double c = p1 / q1;
138         int n = 0;
139         double relativeError = Double.MAX_VALUE;
140         while (n < maxIterations && relativeError > epsilon) {
141             ++n;
142             double a = getA(n, x);
143             double b = getB(n, x);
144             double p2 = a * p1 + b * p0;
145             double q2 = a * q1 + b * q0;
146             if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
147                 // need to scale
148                 if (a != 0.0) {
149                     p2 = p1 + (b / a * p0);
150                     q2 = q1 + (b / a * q0);
151                 } else if (b != 0) {
152                     p2 = (a / b * p1) + p0;
153                     q2 = (a / b * q1) + q0;
154                 } else {
155                     // can not scale an convergent is unbounded.
156                     throw new ConvergenceException(
157                         "Continued fraction convergents diverged to +/- infinity for value {0}",
158                         new Object[] { new Double(x) });
159                 }
160             }
161             double r = p2 / q2;
162             relativeError = Math.abs(r / c - 1.0);
163                 
164             // prepare for next iteration
165             c = p2 / q2;
166             p0 = p1;
167             p1 = p2;
168             q0 = q1;
169             q1 = q2;
170         }
171 
172         if (n >= maxIterations) {
173             throw new MaxIterationsExceededException(maxIterations,
174                 "Continued fraction convergents failed to converge for value {0}",
175                 new Object[] { new Double(x) });
176         }
177 
178         return c;
179     }
180 }