View Javadoc

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