net.openhft.chronicle.hash.impl.util.math.ContinuedFraction Maven / Gradle / Ivy

Go to download
/*
 * Copyright 2012-2018 Chronicle Map Contributors
 * Copyright 2001-2015 The Apache Software Foundation
 * Copyright 2010-2012 CS Systèmes d'Information
 *
 * Licensed 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 net.openhft.chronicle.hash.impl.util.math;

/**
 * Provides a generic means to evaluate continued fractions.  Subclasses simply
 * provided the a and b coefficients to evaluate the continued fraction.
 * 
 * 

 * References:
 * 

 * 
 * Continued Fraction
 * 
 * 
 */
abstract class ContinuedFraction {

    /**
     * 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.
     * 
     * The implementation of this method is based on the modified Lentz algorithm as described
     * on page 18 ff. in:
     * 

     * 
     * I. J. Thompson,  A. R. Barnett. "Coulomb and Bessel Functions of Complex Arguments and Order."
     * 
     * http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf
     * 
     * 
     * Note: the implementation uses the terms a_i and b_i as defined in
     * Continued Fraction @ MathWorld.
     * 
     *
     * @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 IllegalStateException if the algorithm fails to converge.
     * @throws IllegalStateException if maximal number of iterations is reached
     */
    public double evaluate(double x, double epsilon, int maxIterations) {
        final double small = 1e-50;
        double hPrev = getA(0, x);

        // use the value of small as epsilon criteria for zero checks
        if (Precision.isEquals(hPrev, 0.0, small)) {
            hPrev = small;
        }

        int n = 1;
        double dPrev = 0.0;
        double cPrev = hPrev;
        double hN = hPrev;

        while (n < maxIterations) {
            final double a = getA(n, x);
            final double b = getB(n, x);

            double dN = a + b * dPrev;
            if (Precision.isEquals(dN, 0.0, small)) {
                dN = small;
            }
            double cN = a + b / cPrev;
            if (Precision.isEquals(cN, 0.0, small)) {
                cN = small;
            }

            dN = 1 / dN;
            final double deltaN = cN * dN;
            hN = hPrev * deltaN;

            if (Double.isInfinite(hN)) {
                throw new IllegalStateException(
                        "Continued fraction convergents diverged to +/- infinity for value " + x);
            }
            if (Double.isNaN(hN)) {
                throw new IllegalStateException(
                        "Continued fraction diverged to NaN for value " + x);
            }

            if (Math.abs(deltaN - 1.0) < epsilon) {
                break;
            }

            dPrev = dN;
            cPrev = cN;
            hPrev = hN;
            n++;
        }

        if (n >= maxIterations) {
            throw new IllegalStateException(
                    "Continued fraction convergents failed to converge (in less than " +
                            maxIterations + " iterations) for value " + x);
        }

        return hN;
    }

}