/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.util;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* Provides a generic means to evaluate continued fractions. Subclasses simply
* provided the a and b coefficients to evaluate the continued fraction.
*
*
* References:
*
*
*
* @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
*/
public abstract class ContinuedFraction {
/** Maximum allowed numerical error. */
private static final double DEFAULT_EPSILON = 10e-9;
/**
* Default constructor.
*/
protected ContinuedFraction() {
super();
}
/**
* Access the n-th a coefficient of the continued fraction. Since a can be
* a function of the evaluation point, x, that is passed in as well.
* @param n the coefficient index to retrieve.
* @param x the evaluation point.
* @return the n-th a coefficient.
*/
protected abstract double getA(int n, double x);
/**
* Access the n-th b coefficient of the continued fraction. Since b can be
* a function of the evaluation point, x, that is passed in as well.
* @param n the coefficient index to retrieve.
* @param x the evaluation point.
* @return the n-th b coefficient.
*/
protected abstract double getB(int n, double x);
/**
* Evaluates the continued fraction at the value x.
* @param x the evaluation point.
* @return the value of the continued fraction evaluated at x.
* @throws MathException if the algorithm fails to converge.
*/
public double evaluate(double x) throws MathException {
return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
}
/**
* Evaluates the continued fraction at the value x.
* @param x the evaluation point.
* @param epsilon maximum error allowed.
* @return the value of the continued fraction evaluated at x.
* @throws MathException if the algorithm fails to converge.
*/
public double evaluate(double x, double epsilon) throws MathException {
return evaluate(x, epsilon, Integer.MAX_VALUE);
}
/**
* Evaluates the continued fraction at the value x.
* @param x the evaluation point.
* @param maxIterations maximum number of convergents
* @return the value of the continued fraction evaluated at x.
* @throws MathException if the algorithm fails to converge.
*/
public double evaluate(double x, int maxIterations) throws MathException {
return evaluate(x, DEFAULT_EPSILON, maxIterations);
}
/**
*
* Evaluates the continued fraction at the value x.
*
*
*
* The implementation of this method is based on equations 14-17 of:
*
* The recurrence relationship defined in those equations can result in
* very large intermediate results which can result in numerical overflow.
* As a means to combat these overflow conditions, the intermediate results
* are scaled whenever they threaten to become numerically unstable.
*
* @param x the evaluation point.
* @param epsilon maximum error allowed.
* @param maxIterations maximum number of convergents
* @return the value of the continued fraction evaluated at x.
* @throws MathException if the algorithm fails to converge.
*/
public double evaluate(double x, double epsilon, int maxIterations)
throws MathException
{
double p0 = 1.0;
double p1 = getA(0, x);
double q0 = 0.0;
double q1 = 1.0;
double c = p1 / q1;
int n = 0;
double relativeError = Double.MAX_VALUE;
while (n < maxIterations && relativeError > epsilon) {
++n;
double a = getA(n, x);
double b = getB(n, x);
double p2 = a * p1 + b * p0;
double q2 = a * q1 + b * q0;
boolean infinite = false;
if (Double.isInfinite(p2) || Double.isInfinite(q2)) {
/*
* Need to scale. Try successive powers of the larger of a or b
* up to 5th power. Throw ConvergenceException if one or both
* of p2, q2 still overflow.
*/
double scaleFactor = 1d;
double lastScaleFactor = 1d;
final int maxPower = 5;
final double scale = FastMath.max(a,b);
if (scale <= 0) { // Can't scale
throw new ConvergenceException(
LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
x);
}
infinite = true;
for (int i = 0; i < maxPower; i++) {
lastScaleFactor = scaleFactor;
scaleFactor *= scale;
if (a != 0.0 && a > b) {
p2 = p1 / lastScaleFactor + (b / scaleFactor * p0);
q2 = q1 / lastScaleFactor + (b / scaleFactor * q0);
} else if (b != 0) {
p2 = (a / scaleFactor * p1) + p0 / lastScaleFactor;
q2 = (a / scaleFactor * q1) + q0 / lastScaleFactor;
}
infinite = Double.isInfinite(p2) || Double.isInfinite(q2);
if (!infinite) {
break;
}
}
}
if (infinite) {
// Scaling failed
throw new ConvergenceException(
LocalizedFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE,
x);
}
double r = p2 / q2;
if (Double.isNaN(r)) {
throw new ConvergenceException(
LocalizedFormats.CONTINUED_FRACTION_NAN_DIVERGENCE,
x);
}
relativeError = FastMath.abs(r / c - 1.0);
// prepare for next iteration
c = p2 / q2;
p0 = p1;
p1 = p2;
q0 = q1;
q1 = q2;
}
if (n >= maxIterations) {
throw new MaxIterationsExceededException(maxIterations,
LocalizedFormats.NON_CONVERGENT_CONTINUED_FRACTION,
x);
}
return c;
}
}