PolynomialFunctionNewtonForm.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.analysis.polynomials;

  18. import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure;
  19. import org.apache.commons.math4.legacy.analysis.differentiation.UnivariateDifferentiableFunction;
  20. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  21. import org.apache.commons.math4.legacy.exception.NoDataException;
  22. import org.apache.commons.math4.legacy.exception.NullArgumentException;
  23. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;

  24. /**
  25.  * Implements the representation of a real polynomial function in
  26.  * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>,
  27.  * ISBN 0070124477, chapter 2.
  28.  * <p>
  29.  * The formula of polynomial in Newton form is
  30.  *     p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
  31.  *            a[n](x-c[0])(x-c[1])...(x-c[n-1])
  32.  * Note that the length of a[] is one more than the length of c[]</p>
  33.  *
  34.  * @since 1.2
  35.  */
  36. public class PolynomialFunctionNewtonForm implements UnivariateDifferentiableFunction {

  37.     /**
  38.      * The coefficients of the polynomial, ordered by degree -- i.e.
  39.      * coefficients[0] is the constant term and coefficients[n] is the
  40.      * coefficient of x^n where n is the degree of the polynomial.
  41.      */
  42.     private double[] coefficients;

  43.     /**
  44.      * Centers of the Newton polynomial.
  45.      */
  46.     private final double[] c;

  47.     /**
  48.      * When all c[i] = 0, a[] becomes normal polynomial coefficients,
  49.      * i.e. a[i] = coefficients[i].
  50.      */
  51.     private final double[] a;

  52.     /**
  53.      * Whether the polynomial coefficients are available.
  54.      */
  55.     private boolean coefficientsComputed;

  56.     /**
  57.      * Construct a Newton polynomial with the given a[] and c[]. The order of
  58.      * centers are important in that if c[] shuffle, then values of a[] would
  59.      * completely change, not just a permutation of old a[].
  60.      * <p>
  61.      * The constructor makes copy of the input arrays and assigns them.</p>
  62.      *
  63.      * @param a Coefficients in Newton form formula.
  64.      * @param c Centers.
  65.      * @throws NullArgumentException if any argument is {@code null}.
  66.      * @throws NoDataException if any array has zero length.
  67.      * @throws DimensionMismatchException if the size difference between
  68.      * {@code a} and {@code c} is not equal to 1.
  69.      */
  70.     public PolynomialFunctionNewtonForm(double[] a, double[] c)
  71.         throws NullArgumentException, NoDataException, DimensionMismatchException {

  72.         verifyInputArray(a, c);
  73.         this.a = new double[a.length];
  74.         this.c = new double[c.length];
  75.         System.arraycopy(a, 0, this.a, 0, a.length);
  76.         System.arraycopy(c, 0, this.c, 0, c.length);
  77.         coefficientsComputed = false;
  78.     }

  79.     /**
  80.      * Calculate the function value at the given point.
  81.      *
  82.      * @param z Point at which the function value is to be computed.
  83.      * @return the function value.
  84.      */
  85.     @Override
  86.     public double value(double z) {
  87.        return evaluate(a, c, z);
  88.     }

  89.     /**
  90.      * {@inheritDoc}
  91.      * @since 3.1
  92.      */
  93.     @Override
  94.     public DerivativeStructure value(final DerivativeStructure t) {
  95.         verifyInputArray(a, c);

  96.         final int n = c.length;
  97.         DerivativeStructure value = new DerivativeStructure(t.getFreeParameters(), t.getOrder(), a[n]);
  98.         for (int i = n - 1; i >= 0; i--) {
  99.             value = t.subtract(c[i]).multiply(value).add(a[i]);
  100.         }

  101.         return value;
  102.     }

  103.     /**
  104.      * Returns the degree of the polynomial.
  105.      *
  106.      * @return the degree of the polynomial
  107.      */
  108.     public int degree() {
  109.         return c.length;
  110.     }

  111.     /**
  112.      * Returns a copy of coefficients in Newton form formula.
  113.      * <p>
  114.      * Changes made to the returned copy will not affect the polynomial.</p>
  115.      *
  116.      * @return a fresh copy of coefficients in Newton form formula
  117.      */
  118.     public double[] getNewtonCoefficients() {
  119.         double[] out = new double[a.length];
  120.         System.arraycopy(a, 0, out, 0, a.length);
  121.         return out;
  122.     }

  123.     /**
  124.      * Returns a copy of the centers array.
  125.      * <p>
  126.      * Changes made to the returned copy will not affect the polynomial.</p>
  127.      *
  128.      * @return a fresh copy of the centers array.
  129.      */
  130.     public double[] getCenters() {
  131.         double[] out = new double[c.length];
  132.         System.arraycopy(c, 0, out, 0, c.length);
  133.         return out;
  134.     }

  135.     /**
  136.      * Returns a copy of the coefficients array.
  137.      * <p>
  138.      * Changes made to the returned copy will not affect the polynomial.</p>
  139.      *
  140.      * @return a fresh copy of the coefficients array.
  141.      */
  142.     public double[] getCoefficients() {
  143.         if (!coefficientsComputed) {
  144.             computeCoefficients();
  145.         }
  146.         double[] out = new double[coefficients.length];
  147.         System.arraycopy(coefficients, 0, out, 0, coefficients.length);
  148.         return out;
  149.     }

  150.     /**
  151.      * Evaluate the Newton polynomial using nested multiplication. It is
  152.      * also called <a href="http://mathworld.wolfram.com/HornersRule.html">
  153.      * Horner's Rule</a> and takes O(N) time.
  154.      *
  155.      * @param a Coefficients in Newton form formula.
  156.      * @param c Centers.
  157.      * @param z Point at which the function value is to be computed.
  158.      * @return the function value.
  159.      * @throws NullArgumentException if any argument is {@code null}.
  160.      * @throws NoDataException if any array has zero length.
  161.      * @throws DimensionMismatchException if the size difference between
  162.      * {@code a} and {@code c} is not equal to 1.
  163.      */
  164.     public static double evaluate(double[] a, double[] c, double z)
  165.         throws NullArgumentException, DimensionMismatchException, NoDataException {
  166.         verifyInputArray(a, c);

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

  172.         return value;
  173.     }

  174.     /**
  175.      * Calculate the normal polynomial coefficients given the Newton form.
  176.      * It also uses nested multiplication but takes O(N^2) time.
  177.      */
  178.     protected void computeCoefficients() {
  179.         final int n = degree();

  180.         coefficients = new double[n+1];
  181.         for (int i = 0; i <= n; i++) {
  182.             coefficients[i] = 0.0;
  183.         }

  184.         coefficients[0] = a[n];
  185.         for (int i = n-1; i >= 0; i--) {
  186.             for (int j = n-i; j > 0; j--) {
  187.                 coefficients[j] = coefficients[j-1] - c[i] * coefficients[j];
  188.             }
  189.             coefficients[0] = a[i] - c[i] * coefficients[0];
  190.         }

  191.         coefficientsComputed = true;
  192.     }

  193.     /**
  194.      * Verifies that the input arrays are valid.
  195.      * <p>
  196.      * The centers must be distinct for interpolation purposes, but not
  197.      * for general use. Thus it is not verified here.</p>
  198.      *
  199.      * @param a the coefficients in Newton form formula
  200.      * @param c the centers
  201.      * @throws NullArgumentException if any argument is {@code null}.
  202.      * @throws NoDataException if any array has zero length.
  203.      * @throws DimensionMismatchException if the size difference between
  204.      * {@code a} and {@code c} is not equal to 1.
  205.      * @see org.apache.commons.math4.legacy.analysis.interpolation.DividedDifferenceInterpolator#computeDividedDifference(double[],
  206.      * double[])
  207.      */
  208.     protected static void verifyInputArray(double[] a, double[] c)
  209.         throws NullArgumentException, NoDataException, DimensionMismatchException {
  210.         NullArgumentException.check(a);
  211.         NullArgumentException.check(c);
  212.         if (a.length == 0 || c.length == 0) {
  213.             throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
  214.         }
  215.         if (a.length != c.length + 1) {
  216.             throw new DimensionMismatchException(LocalizedFormats.ARRAY_SIZES_SHOULD_HAVE_DIFFERENCE_1,
  217.                                                  a.length, c.length);
  218.         }
  219.     }
  220. }