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: 1.14 $ $Date: 2004/06/23 16:26:16 $
36   */
37  public abstract class ContinuedFraction implements Serializable {
38      
39      /*** Serialization UID */
40      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      * Evaluates the continued fraction at the value x.
104      * 
105      * The implementation of this method is based on:
106      * <ul>
107      * <li>O. E-gecio-glu, C . K. Koc, J. Rifa i Coma,
108      * <a href="http://citeseer.ist.psu.edu/egecioglu91fast.html">
109      * On Fast Computation of Continued Fractions</a>, Computers Math. Applic.,
110      * 21(2--3), 1991, 167--169.</li>
111      * </ul>
112      * 
113      * @param x the evaluation point.
114      * @param epsilon maximum error allowed.
115      * @param maxIterations maximum number of convergents
116      * @return the value of the continued fraction evaluated at x. 
117      * @throws MathException if the algorithm fails to converge.
118      */
119     public double evaluate(double x, double epsilon, int maxIterations)
120         throws MathException
121     {
122         double[][] f = new double[2][2];
123         double[][] a = new double[2][2];
124         double[][] an = new double[2][2];
125 
126         a[0][0] = getA(0, x);
127         a[0][1] = 1.0;
128         a[1][0] = 1.0;
129         a[1][1] = 0.0;
130 
131         return evaluate(1, x, a, an, f, epsilon, maxIterations);
132     }
133 
134     /***
135      * Evaluates the n-th convergent, fn = pn / qn, for this continued fraction
136      * at the value x.
137      * @param n the convergent to compute.
138      * @param x the evaluation point.
139      * @param a (n-1)-th convergent matrix.  (Input)
140      * @param an the n-th coefficient matrix. (Output)
141      * @param f the n-th convergent matrix. (Output)
142      * @param epsilon maximum error allowed.
143      * @param maxIterations maximum number of convergents
144      * @return the value of the the n-th convergent for this continued fraction
145      *         evaluated at x. 
146      * @throws MathException if the algorithm fails to converge.
147      */
148     private double evaluate(
149         int n,
150         double x,
151         double[][] a,
152         double[][] an,
153         double[][] f,
154         double epsilon,
155         int maxIterations) throws MathException 
156     {
157         double ret;
158 
159         // create next matrix
160         an[0][0] = getA(n, x);
161         an[0][1] = 1.0;
162         an[1][0] = getB(n, x);
163         an[1][1] = 0.0;
164 
165         // multiply a and an, save as f
166         f[0][0] = (a[0][0] * an[0][0]) + (a[0][1] * an[1][0]);
167         f[0][1] = (a[0][0] * an[0][1]) + (a[0][1] * an[1][1]);
168         f[1][0] = (a[1][0] * an[0][0]) + (a[1][1] * an[1][0]);
169         f[1][1] = (a[1][0] * an[0][1]) + (a[1][1] * an[1][1]);
170 
171         // determine if we're close enough
172         if (Math.abs((f[0][0] * f[1][1]) - (f[1][0] * f[0][1])) <
173             Math.abs(epsilon * f[1][0] * f[1][1]))
174         {
175             ret = f[0][0] / f[1][0];
176         } else {
177             if (n >= maxIterations) {
178                 throw new ConvergenceException(
179                     "Continued fraction convergents failed to converge.");
180             }
181             // compute next
182             ret = evaluate(n + 1, x, f /* new a */
183             , an /* reuse an */
184             , a /* new f */
185             , epsilon, maxIterations);
186         }
187 
188         return ret;
189     }
190 }