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 org.apache.commons.math.ConvergenceException;
20  import org.apache.commons.math.MathException;
21  import org.apache.commons.math.MaxIterationsExceededException;
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: 778522 $ $Date: 2009-05-25 18:24:50 -0400 (Mon, 25 May 2009) $
36   */
37  public abstract class ContinuedFraction {
38      
39      /** Maximum allowed numerical error. */
40      private static final double DEFAULT_EPSILON = 10e-9;
41  
42      /**
43       * Default constructor.
44       */
45      protected ContinuedFraction() {
46          super();
47      }
48  
49      /**
50       * Access the n-th a coefficient of the continued fraction.  Since a can be
51       * a function of the evaluation point, x, that is passed in as well.
52       * @param n the coefficient index to retrieve.
53       * @param x the evaluation point.
54       * @return the n-th a coefficient.
55       */
56      protected abstract double getA(int n, double x);
57  
58      /**
59       * Access the n-th b coefficient of the continued fraction.  Since b can be
60       * a function of the evaluation point, x, that is passed in as well.
61       * @param n the coefficient index to retrieve.
62       * @param x the evaluation point.
63       * @return the n-th b coefficient.
64       */
65      protected abstract double getB(int n, double x);
66  
67      /**
68       * Evaluates the continued fraction at the value x.
69       * @param x the evaluation point.
70       * @return the value of the continued fraction evaluated at x. 
71       * @throws MathException if the algorithm fails to converge.
72       */
73      public double evaluate(double x) throws MathException {
74          return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
75      }
76  
77      /**
78       * Evaluates the continued fraction at the value x.
79       * @param x the evaluation point.
80       * @param epsilon maximum error allowed.
81       * @return the value of the continued fraction evaluated at x. 
82       * @throws MathException if the algorithm fails to converge.
83       */
84      public double evaluate(double x, double epsilon) throws MathException {
85          return evaluate(x, epsilon, Integer.MAX_VALUE);
86      }
87  
88      /**
89       * Evaluates the continued fraction at the value x.
90       * @param x the evaluation point.
91       * @param maxIterations maximum number of convergents
92       * @return the value of the continued fraction evaluated at x. 
93       * @throws MathException if the algorithm fails to converge.
94       */
95      public double evaluate(double x, int maxIterations) throws MathException {
96          return evaluate(x, DEFAULT_EPSILON, maxIterations);
97      }
98  
99      /**
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 }