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 }