FieldContinuedFraction

/*
 * Licensed to the Hipparchus project under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The Hipparchus project 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
 *
 *      https://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.hipparchus.util;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalStateException;

/**
 * Provides a generic means to evaluate continued fractions.  Subclasses simply
 * provided the a and b coefficients to evaluate the continued fraction.
 * <p>
 * References:
 * <ul>
 * <li><a href="http://mathworld.wolfram.com/ContinuedFraction.html">
 * Continued Fraction</a></li>
 * </ul>
 *
 */
public abstract class FieldContinuedFraction {
    /** Maximum allowed numerical error. */
    private static final double DEFAULT_EPSILON = 10e-9;

    /**
     * Default constructor.
     */
    protected FieldContinuedFraction() {
        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.
     * @param <T> type of the field elements.
     * @return the n-th a coefficient.
     */
    public abstract <T extends CalculusFieldElement<T>> T getA(int n, T 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.
     * @param <T> type of the field elements.
     * @return the n-th b coefficient.
     */
    public abstract <T extends CalculusFieldElement<T>> T getB(int n, T x);

    /**
     * Evaluates the continued fraction at the value x.
     * @param x the evaluation point.
     * @param <T> type of the field elements.
     * @return the value of the continued fraction evaluated at x.
     * @throws MathIllegalStateException if the algorithm fails to converge.
     */
    public <T extends CalculusFieldElement<T>> T evaluate(T x) throws MathIllegalStateException {
        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.
     * @param <T> type of the field elements.
     * @return the value of the continued fraction evaluated at x.
     * @throws MathIllegalStateException if the algorithm fails to converge.
     */
    public <T extends CalculusFieldElement<T>> T evaluate(T x, double epsilon) throws MathIllegalStateException {
        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
     * @param <T> type of the field elements.
     * @return the value of the continued fraction evaluated at x.
     * @throws MathIllegalStateException if the algorithm fails to converge.
     * @throws MathIllegalStateException if maximal number of iterations is reached
     */
    public <T extends CalculusFieldElement<T>> T evaluate(T x, int maxIterations)
        throws MathIllegalStateException {
        return evaluate(x, DEFAULT_EPSILON, maxIterations);
    }

    /**
     * Evaluates the continued fraction at the value x.
     * <p>
     * The implementation of this method is based on the modified Lentz algorithm as described
     * on page 18 ff. in:
     * </p>
     * <ul>
     *   <li>
     *   I. J. Thompson,  A. R. Barnett. "Coulomb and Bessel Functions of Complex Arguments and Order."
     *   <a target="_blank" href="http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf">
     *   http://www.fresco.org.uk/papers/Thompson-JCP64p490.pdf</a>
     *   </li>
     * </ul>
     * <p>
     * <b>Note:</b> the implementation uses the terms a<sub>i</sub> and b<sub>i</sub> as defined in
     * <a href="http://mathworld.wolfram.com/ContinuedFraction.html">Continued Fraction @ MathWorld</a>.
     * </p>
     *
     * @param x the evaluation point.
     * @param epsilon maximum error allowed.
     * @param maxIterations maximum number of convergents
     * @param <T> type of the field elements.
     * @return the value of the continued fraction evaluated at x.
     * @throws MathIllegalStateException if the algorithm fails to converge.
     * @throws MathIllegalStateException if maximal number of iterations is reached
     */
    public <T extends CalculusFieldElement<T>> T evaluate(T x, double epsilon, int maxIterations)
        throws MathIllegalStateException {
        final T zero = x.getField().getZero();
        final T one  = x.getField().getOne();

        final double small      = 1e-50;
        final T      smallField = one.multiply(small);

        T hPrev = getA(0, x);

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

        int n     = 1;
        T   dPrev = zero;
        T   cPrev = hPrev;
        T   hN    = hPrev;

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

            T dN = a.add(b.multiply(dPrev));
            if (Precision.equals(dN.getReal(), 0.0, small)) {
                dN = smallField;
            }
            T cN = a.add(b.divide(cPrev));
            if (Precision.equals(cN.getReal(), 0.0, small)) {
                cN = smallField;
            }

            dN = dN.reciprocal();
            final T deltaN = cN.multiply(dN);
            hN = hPrev.multiply(deltaN);

            if (hN.isInfinite()) {
                throw new MathIllegalStateException(LocalizedCoreFormats.CONTINUED_FRACTION_INFINITY_DIVERGENCE, x);
            }
            if (hN.isNaN()) {
                throw new MathIllegalStateException(LocalizedCoreFormats.CONTINUED_FRACTION_NAN_DIVERGENCE, x);
            }

            if (deltaN.subtract(1.0).abs().getReal() < epsilon) {
                break;
            }

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

        if (n >= maxIterations) {
            throw new MathIllegalStateException(LocalizedCoreFormats.NON_CONVERGENT_CONTINUED_FRACTION,
                                                maxIterations, x);
        }

        return hN;
    }

}