Coverage Report - org.apache.commons.math.util.ContinuedFraction

Classes in this File Line Coverage Branch Coverage Complexity
ContinuedFraction
81% 
80% 
2

 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$ $Date: 2005-08-22 19:27:17 -0700 (Mon, 22 Aug 2005) $
 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  6282
         super();
 50  6282
     }
 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  0
         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  2
         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  0
         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  6282
             double p0 = 1.0;
 132  6282
             double p1 = getA(0, x);
 133  6282
             double q0 = 0.0;
 134  6282
             double q1 = 1.0;
 135  6282
             double c = p1 / q1;
 136  6282
             int n = 0;
 137  6282
             double relativeError = Double.MAX_VALUE;
 138  65950
             while (n < maxIterations && relativeError > epsilon) {
 139  59668
                     ++n;
 140  59668
                     double a = getA(n, x);
 141  59668
                     double b = getB(n, x);
 142  59668
                           double p2 = a * p1 + b * p0;
 143  59668
                            double q2 = a * q1 + b * q0;
 144  59668
                            if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
 145  
                                    // need to scale
 146  628
                                    if (a != 0.0) {
 147  628
                                            p2 = p1 + (b / a * p0);
 148  628
                                            q2 = q1 + (b / a * q0);
 149  0
                                    } else if (b != 0) {
 150  0
                                            p2 = (a / b * p1) + p0;
 151  0
                                            q2 = (a / b * q1) + q0;
 152  
                                    } else {
 153  
                                            // can not scale an convergent is unbounded.
 154  0
                                throw new ConvergenceException(
 155  
                                    "Continued fraction convergents diverged to +/- " +
 156  
                                    "infinity.");
 157  
                                    }
 158  
                            }
 159  59668
                            double r = p2 / q2;
 160  59668
                            relativeError = Math.abs(r / c - 1.0);
 161  
                             
 162  
                            // prepare for next iteration
 163  59668
                            c = p2 / q2;
 164  59668
                            p0 = p1;
 165  59668
                            p1 = p2;
 166  59668
                            q0 = q1;
 167  59668
                            q1 = q2;
 168  
             }
 169  
 
 170  6282
             if (n >= maxIterations) {
 171  0
             throw new ConvergenceException(
 172  
                 "Continued fraction convergents failed to converge.");
 173  
         }
 174  
 
 175  6282
         return c;
 176  
     }
 177  
 }