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 }