FieldPolynomialFunctionLagrangeForm.java

/*
 * 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.analysis.polynomials;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;

/**
 * Implements the representation of a real polynomial function in
 * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
 * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
 * Analysis</b>, ISBN 038795452X, chapter 2.
 * <p>
 * The approximated function should be smooth enough for Lagrange polynomial
 * to work well. Otherwise, consider using splines instead.</p>
 * @see PolynomialFunctionLagrangeForm
 * @since 4.0
 * @param <T> type of the field elements
 */
public class FieldPolynomialFunctionLagrangeForm<T extends CalculusFieldElement<T>>
        implements CalculusFieldUnivariateFunction<T> {
    /**
     * The coefficients of the polynomial, ordered by degree -- i.e.
     * coefficients[0] is the constant term and coefficients[n] is the
     * coefficient of x^n where n is the degree of the polynomial.
     */
    private T[] coefficients;
    /**
     * Interpolating points (abscissas).
     */
    private final T[] x;
    /**
     * Function values at interpolating points.
     */
    private final T[] y;
    /**
     * Whether the polynomial coefficients are available.
     */
    private boolean coefficientsComputed;

    /**
     * Construct a Lagrange polynomial with the given abscissas and function
     * values. The order of interpolating points is important.
     * <p>
     * The constructor makes copy of the input arrays and assigns them.</p>
     *
     * @param x interpolating points
     * @param y function values at interpolating points
     * @throws MathIllegalArgumentException if the array lengths are different.
     * @throws MathIllegalArgumentException if the number of points is less than 2.
     * @throws MathIllegalArgumentException if two abscissae have the same value.
     * @throws MathIllegalArgumentException if the abscissae are not sorted.
     */
    public FieldPolynomialFunctionLagrangeForm(final T[] x, final T[] y)
        throws MathIllegalArgumentException {
        this.x = x.clone();
        this.y = y.clone();
        coefficientsComputed = false;

        MathArrays.checkEqualLength(x, y);
        if (x.length < 2) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
        }
        MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, true);
    }

    /**
     * Calculate the function value at the given point.
     *
     * @param z Point at which the function value is to be computed.
     * @return the function value.
     * @throws MathIllegalArgumentException if {@code x} and {@code y} have
     * different lengths.
     * @throws MathIllegalArgumentException
     * if {@code x} is not sorted in strictly increasing order.
     * @throws MathIllegalArgumentException if the size of {@code x} is less
     * than 2.
     */
    @Override
    public T value(final T z) {
        int nearest = 0;
        final int n = x.length;
        final T[] c = y.clone();
        final T[] d = c.clone();
        double minDist = Double.POSITIVE_INFINITY;
        for (int i = 0; i < n; i++) {
            // find out the abscissa closest to z
            final double dist = FastMath.abs(z.subtract(x[i])).getReal();
            if (dist < minDist) {
                nearest = i;
                minDist = dist;
            }
        }

        // initial approximation to the function value at z
        T value = y[nearest];

        for (int i = 1; i < n; i++) {
            for (int j = 0; j < n-i; j++) {
                final T tc = x[j].subtract(z);
                final T td = x[i+j].subtract(z);
                final T divider = x[j].subtract(x[i+j]);
                // update the difference arrays
                final T w = (c[j+1].subtract(d[j])).divide(divider);
                c[j] = tc.multiply(w);
                d[j] = td.multiply(w);
            }
            // sum up the difference terms to get the final value
            if (nearest < 0.5*(n-i+1)) {
                value = value.add(c[nearest]);    // fork down
            } else {
                nearest--;
                value = value.add(d[nearest]);    // fork up
            }
        }

        return value;
    }

    /**
     * Returns the degree of the polynomial.
     *
     * @return the degree of the polynomial
     */
    public int degree() {
        return x.length - 1;
    }

    /**
     * Returns a copy of the interpolating points array.
     * <p>
     * Changes made to the returned copy will not affect the polynomial.</p>
     *
     * @return a fresh copy of the interpolating points array
     */
    public T[] getInterpolatingPoints() {
        return x.clone();
    }

    /**
     * Returns a copy of the interpolating values array.
     * <p>
     * Changes made to the returned copy will not affect the polynomial.</p>
     *
     * @return a fresh copy of the interpolating values array
     */
    public T[] getInterpolatingValues() {
        return y.clone();
    }

    /**
     * Returns a copy of the coefficients array.
     * <p>
     * Changes made to the returned copy will not affect the polynomial.</p>
     * <p>
     * Note that coefficients computation can be ill-conditioned. Use with caution
     * and only when it is necessary.</p>
     *
     * @return a fresh copy of the coefficients array
     */
    public T[] getCoefficients() {
        if (!coefficientsComputed) {
            computeCoefficients();
        }
        return coefficients.clone();
    }

    /**
     * Calculate the coefficients of Lagrange polynomial from the
     * interpolation data. It takes O(n^2) time.
     * Note that this computation can be ill-conditioned: Use with caution
     * and only when it is necessary.
     */
    protected void computeCoefficients() {
        final int n = degree() + 1;
        final Field<T> field = x[0].getField();
        coefficients = MathArrays.buildArray(field, n);

        // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
        final T[] c = MathArrays.buildArray(field, n + 1);
        c[0] = field.getOne();
        for (int i = 0; i < n; i++) {
            for (int j = i; j > 0; j--) {
                c[j] = c[j-1].subtract(c[j].multiply(x[i]));
            }
            c[0] = c[0].multiply(x[i].negate());
            c[i+1] = field.getOne();
        }

        final T[] tc = MathArrays.buildArray(field, n);
        for (int i = 0; i < n; i++) {
            // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
            T d = field.getOne();
            for (int j = 0; j < n; j++) {
                if (i != j) {
                    d = d.multiply(x[i].subtract(x[j]));
                }
            }
            final T t = y[i].divide(d);
            // Lagrange polynomial is the sum of n terms, each of which is a
            // polynomial of degree n-1. tc[] are the coefficients of the i-th
            // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
            tc[n-1] = c[n];     // actually c[n] = 1
            coefficients[n-1] = coefficients[n-1].add(t.multiply(tc[n-1]));
            for (int j = n-2; j >= 0; j--) {
                tc[j] = c[j+1].add(tc[j+1].multiply(x[i]));
                coefficients[j] = coefficients[j].add(t.multiply(tc[j]));
            }
        }

        coefficientsComputed = true;
    }
}