1 files changed, 245 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionNewtonForm.java b/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionNewtonForm.java
new file mode 100644
index 0000000..fc2f1fd
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionNewtonForm.java
@@ -0,0 +1,245 @@
+/*
+ * 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.math3.analysis.polynomials;
+
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
+import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.NoDataException;
+import org.apache.commons.math3.exception.NullArgumentException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.util.MathUtils;
+
+/**
+ * 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>
+ *
+ * @since 1.2
+ */
+public class PolynomialFunctionNewtonForm implements UnivariateDifferentiableFunction {
+
+    /**
+     * 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[];
+
+    /**
+     * Centers of the Newton polynomial.
+     */
+    private final double c[];
+
+    /**
+     * When all c[i] = 0, a[] becomes normal polynomial coefficients,
+     * i.e. a[i] = coefficients[i].
+     */
+    private final double a[];
+
+    /**
+     * 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 Coefficients in Newton form formula.
+     * @param c Centers.
+     * @throws NullArgumentException if any argument is {@code null}.
+     * @throws NoDataException if any array has zero length.
+     * @throws DimensionMismatchException if the size difference between
+     * {@code a} and {@code c} is not equal to 1.
+     */
+    public PolynomialFunctionNewtonForm(double a[], double c[])
+        throws NullArgumentException, NoDataException, DimensionMismatchException {
+
+        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 Point at which the function value is to be computed.
+     * @return the function value.
+     */
+    public double value(double z) {
+       return evaluate(a, c, z);
+    }
+
+    /**
+     * {@inheritDoc}
+     * @since 3.1
+     */
+    public DerivativeStructure value(final DerivativeStructure t) {
+        verifyInputArray(a, c);
+
+        final int n = c.length;
+        DerivativeStructure value = new DerivativeStructure(t.getFreeParameters(), t.getOrder(), a[n]);
+        for (int i = n - 1; i >= 0; i--) {
+            value = t.subtract(c[i]).multiply(value).add(a[i]);
+        }
+
+        return value;
+
+    }
+
+    /**
+     * 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 Coefficients in Newton form formula.
+     * @param c Centers.
+     * @param z Point at which the function value is to be computed.
+     * @return the function value.
+     * @throws NullArgumentException if any argument is {@code null}.
+     * @throws NoDataException if any array has zero length.
+     * @throws DimensionMismatchException if the size difference between
+     * {@code a} and {@code c} is not equal to 1.
+     */
+    public static double evaluate(double a[], double c[], double z)
+        throws NullArgumentException, DimensionMismatchException, NoDataException {
+        verifyInputArray(a, c);
+
+        final 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() {
+        final int n = degree();
+
+        coefficients = new double[n+1];
+        for (int i = 0; i <= n; i++) {
+            coefficients[i] = 0.0;
+        }
+
+        coefficients[0] = a[n];
+        for (int i = n-1; i >= 0; i--) {
+            for (int 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 NullArgumentException if any argument is {@code null}.
+     * @throws NoDataException if any array has zero length.
+     * @throws DimensionMismatchException if the size difference between
+     * {@code a} and {@code c} is not equal to 1.
+     * @see org.apache.commons.math3.analysis.interpolation.DividedDifferenceInterpolator#computeDividedDifference(double[],
+     * double[])
+     */
+    protected static void verifyInputArray(double a[], double c[])
+        throws NullArgumentException, NoDataException, DimensionMismatchException {
+        MathUtils.checkNotNull(a);
+        MathUtils.checkNotNull(c);
+        if (a.length == 0 || c.length == 0) {
+            throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
+        }
+        if (a.length != c.length + 1) {
+            throw new DimensionMismatchException(LocalizedFormats.ARRAY_SIZES_SHOULD_HAVE_DIFFERENCE_1,
+                                                 a.length, c.length);
+        }
+    }
+
+}