← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/util/ContinuedFraction.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/util/ContinuedFraction.java
 
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
}