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

import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.interpolation.DividedDifferenceInterpolator;

/**
 * Implements the representation of a real polynomial function in Newton Form.
 * For reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
 * chapter 2.
 * <p>
 * The formula of polynomial in Newton form is p(x) = a[0] + a[1](x-c[0]) +
 * a[2](x-c[0])(x-c[1]) + ... + a[n](x-c[0])(x-c[1])...(x-c[n-1]) Note that the
 * length of a[] is one more than the length of c[]
 * </p>
 *
 * @version $Revision$ $Date$
 * @since 1.2
 */
public class PolynomialFunctionNewtonForm implements UnivariateRealFunction {

	/**
	 * 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 double coefficients[];

	/**
	 * Members of c[] are called centers of the Newton polynomial. When all c[i]
	 * = 0, a[] becomes normal polynomial coefficients, i.e. a[i] =
	 * coefficients[i].
	 */
	private double a[], c[];

	/**
	 * Whether the polynomial coefficients are available.
	 */
	private boolean coefficientsComputed;

	/**
	 * Construct a Newton polynomial with the given a[] and c[]. The order of
	 * centers are important in that if c[] shuffle, then values of a[] would
	 * completely change, not just a permutation of old a[].
	 * <p>
	 * The constructor makes copy of the input arrays and assigns them.
	 * </p>
	 * 
	 * @param a
	 *            the coefficients in Newton form formula
	 * @param c
	 *            the centers
	 * @throws IllegalArgumentException
	 *             if input arrays are not valid
	 */
	public PolynomialFunctionNewtonForm(double a[], double c[]) throws IllegalArgumentException {

		verifyInputArray(a, c);
		this.a = new double[a.length];
		this.c = new double[c.length];
		System.arraycopy(a, 0, this.a, 0, a.length);
		System.arraycopy(c, 0, this.c, 0, c.length);
		coefficientsComputed = false;
	}

	/**
	 * Calculate the function value at the given point.
	 *
	 * @param z
	 *            the point at which the function value is to be computed
	 * @return the function value
	 * @throws FunctionEvaluationException
	 *             if a runtime error occurs
	 * @see UnivariateRealFunction#value(double)
	 */
	public double value(double z) throws FunctionEvaluationException {
		return evaluate(a, c, z);
	}

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

	/**
	 * Returns a copy of coefficients in Newton form formula.
	 * <p>
	 * Changes made to the returned copy will not affect the polynomial.
	 * </p>
	 * 
	 * @return a fresh copy of coefficients in Newton form formula
	 */
	public double[] getNewtonCoefficients() {
		double[] out = new double[a.length];
		System.arraycopy(a, 0, out, 0, a.length);
		return out;
	}

	/**
	 * Returns a copy of the centers array.
	 * <p>
	 * Changes made to the returned copy will not affect the polynomial.
	 * </p>
	 * 
	 * @return a fresh copy of the centers array
	 */
	public double[] getCenters() {
		double[] out = new double[c.length];
		System.arraycopy(c, 0, out, 0, c.length);
		return out;
	}

	/**
	 * Returns a copy of the coefficients array.
	 * <p>
	 * Changes made to the returned copy will not affect the polynomial.
	 * </p>
	 * 
	 * @return a fresh copy of the coefficients array
	 */
	public double[] getCoefficients() {
		if (!coefficientsComputed) {
			computeCoefficients();
		}
		double[] out = new double[coefficients.length];
		System.arraycopy(coefficients, 0, out, 0, coefficients.length);
		return out;
	}

	/**
	 * Evaluate the Newton polynomial using nested multiplication. It is also
	 * called <a href="http://mathworld.wolfram.com/HornersRule.html"> Horner's
	 * Rule</a> and takes O(N) time.
	 *
	 * @param a
	 *            the coefficients in Newton form formula
	 * @param c
	 *            the centers
	 * @param z
	 *            the point at which the function value is to be computed
	 * @return the function value
	 * @throws FunctionEvaluationException
	 *             if a runtime error occurs
	 * @throws IllegalArgumentException
	 *             if inputs are not valid
	 */
	public static double evaluate(double a[], double c[], double z)
			throws FunctionEvaluationException, IllegalArgumentException {

		verifyInputArray(a, c);

		int n = c.length;
		double value = a[n];
		for (int i = n - 1; i >= 0; i--) {
			value = a[i] + (z - c[i]) * value;
		}

		return value;
	}

	/**
	 * Calculate the normal polynomial coefficients given the Newton form. It
	 * also uses nested multiplication but takes O(N^2) time.
	 */
	protected void computeCoefficients() {
		int i, j, n = degree();

		coefficients = new double[n + 1];
		for (i = 0; i <= n; i++) {
			coefficients[i] = 0.0;
		}

		coefficients[0] = a[n];
		for (i = n - 1; i >= 0; i--) {
			for (j = n - i; j > 0; j--) {
				coefficients[j] = coefficients[j - 1] - c[i] * coefficients[j];
			}
			coefficients[0] = a[i] - c[i] * coefficients[0];
		}

		coefficientsComputed = true;
	}

	/**
	 * Verifies that the input arrays are valid.
	 * <p>
	 * The centers must be distinct for interpolation purposes, but not for
	 * general use. Thus it is not verified here.
	 * </p>
	 * 
	 * @param a
	 *            the coefficients in Newton form formula
	 * @param c
	 *            the centers
	 * @throws IllegalArgumentException
	 *             if not valid
	 * @see DividedDifferenceInterpolator#computeDividedDifference(double[],
	 *      double[])
	 */
	protected static void verifyInputArray(double a[], double c[]) throws IllegalArgumentException {

		if (a.length < 1 || c.length < 1) {
			throw MathRuntimeException.createIllegalArgumentException("empty polynomials coefficients array");
		}
		if (a.length != c.length + 1) {
			throw MathRuntimeException.createIllegalArgumentException(
					"array sizes should have difference 1 ({0} != {1} + 1)", a.length, c.length);
		}
	}
}
