001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.analysis.polynomials;
018
019import org.apache.commons.math4.analysis.differentiation.DerivativeStructure;
020import org.apache.commons.math4.analysis.differentiation.UnivariateDifferentiableFunction;
021import org.apache.commons.math4.exception.DimensionMismatchException;
022import org.apache.commons.math4.exception.NoDataException;
023import org.apache.commons.math4.exception.NullArgumentException;
024import org.apache.commons.math4.exception.util.LocalizedFormats;
025import org.apache.commons.math4.util.MathUtils;
026
027/**
028 * Implements the representation of a real polynomial function in
029 * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>,
030 * ISBN 0070124477, chapter 2.
031 * <p>
032 * The formula of polynomial in Newton form is
033 *     p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
034 *            a[n](x-c[0])(x-c[1])...(x-c[n-1])
035 * Note that the length of a[] is one more than the length of c[]</p>
036 *
037 * @since 1.2
038 */
039public class PolynomialFunctionNewtonForm implements UnivariateDifferentiableFunction {
040
041    /**
042     * The coefficients of the polynomial, ordered by degree -- i.e.
043     * coefficients[0] is the constant term and coefficients[n] is the
044     * coefficient of x^n where n is the degree of the polynomial.
045     */
046    private double coefficients[];
047
048    /**
049     * Centers of the Newton polynomial.
050     */
051    private final double c[];
052
053    /**
054     * When all c[i] = 0, a[] becomes normal polynomial coefficients,
055     * i.e. a[i] = coefficients[i].
056     */
057    private final double a[];
058
059    /**
060     * Whether the polynomial coefficients are available.
061     */
062    private boolean coefficientsComputed;
063
064    /**
065     * Construct a Newton polynomial with the given a[] and c[]. The order of
066     * centers are important in that if c[] shuffle, then values of a[] would
067     * completely change, not just a permutation of old a[].
068     * <p>
069     * The constructor makes copy of the input arrays and assigns them.</p>
070     *
071     * @param a Coefficients in Newton form formula.
072     * @param c Centers.
073     * @throws NullArgumentException if any argument is {@code null}.
074     * @throws NoDataException if any array has zero length.
075     * @throws DimensionMismatchException if the size difference between
076     * {@code a} and {@code c} is not equal to 1.
077     */
078    public PolynomialFunctionNewtonForm(double a[], double c[])
079        throws NullArgumentException, NoDataException, DimensionMismatchException {
080
081        verifyInputArray(a, c);
082        this.a = new double[a.length];
083        this.c = new double[c.length];
084        System.arraycopy(a, 0, this.a, 0, a.length);
085        System.arraycopy(c, 0, this.c, 0, c.length);
086        coefficientsComputed = false;
087    }
088
089    /**
090     * Calculate the function value at the given point.
091     *
092     * @param z Point at which the function value is to be computed.
093     * @return the function value.
094     */
095    @Override
096    public double value(double z) {
097       return evaluate(a, c, z);
098    }
099
100    /**
101     * {@inheritDoc}
102     * @since 3.1
103     */
104    @Override
105    public DerivativeStructure value(final DerivativeStructure t) {
106        verifyInputArray(a, c);
107
108        final int n = c.length;
109        DerivativeStructure value = new DerivativeStructure(t.getFreeParameters(), t.getOrder(), a[n]);
110        for (int i = n - 1; i >= 0; i--) {
111            value = t.subtract(c[i]).multiply(value).add(a[i]);
112        }
113
114        return value;
115
116    }
117
118    /**
119     * Returns the degree of the polynomial.
120     *
121     * @return the degree of the polynomial
122     */
123    public int degree() {
124        return c.length;
125    }
126
127    /**
128     * Returns a copy of coefficients in Newton form formula.
129     * <p>
130     * Changes made to the returned copy will not affect the polynomial.</p>
131     *
132     * @return a fresh copy of coefficients in Newton form formula
133     */
134    public double[] getNewtonCoefficients() {
135        double[] out = new double[a.length];
136        System.arraycopy(a, 0, out, 0, a.length);
137        return out;
138    }
139
140    /**
141     * Returns a copy of the centers array.
142     * <p>
143     * Changes made to the returned copy will not affect the polynomial.</p>
144     *
145     * @return a fresh copy of the centers array.
146     */
147    public double[] getCenters() {
148        double[] out = new double[c.length];
149        System.arraycopy(c, 0, out, 0, c.length);
150        return out;
151    }
152
153    /**
154     * Returns a copy of the coefficients array.
155     * <p>
156     * Changes made to the returned copy will not affect the polynomial.</p>
157     *
158     * @return a fresh copy of the coefficients array.
159     */
160    public double[] getCoefficients() {
161        if (!coefficientsComputed) {
162            computeCoefficients();
163        }
164        double[] out = new double[coefficients.length];
165        System.arraycopy(coefficients, 0, out, 0, coefficients.length);
166        return out;
167    }
168
169    /**
170     * Evaluate the Newton polynomial using nested multiplication. It is
171     * also called <a href="http://mathworld.wolfram.com/HornersRule.html">
172     * Horner's Rule</a> and takes O(N) time.
173     *
174     * @param a Coefficients in Newton form formula.
175     * @param c Centers.
176     * @param z Point at which the function value is to be computed.
177     * @return the function value.
178     * @throws NullArgumentException if any argument is {@code null}.
179     * @throws NoDataException if any array has zero length.
180     * @throws DimensionMismatchException if the size difference between
181     * {@code a} and {@code c} is not equal to 1.
182     */
183    public static double evaluate(double a[], double c[], double z)
184        throws NullArgumentException, DimensionMismatchException, NoDataException {
185        verifyInputArray(a, c);
186
187        final int n = c.length;
188        double value = a[n];
189        for (int i = n - 1; i >= 0; i--) {
190            value = a[i] + (z - c[i]) * value;
191        }
192
193        return value;
194    }
195
196    /**
197     * Calculate the normal polynomial coefficients given the Newton form.
198     * It also uses nested multiplication but takes O(N^2) time.
199     */
200    protected void computeCoefficients() {
201        final int n = degree();
202
203        coefficients = new double[n+1];
204        for (int i = 0; i <= n; i++) {
205            coefficients[i] = 0.0;
206        }
207
208        coefficients[0] = a[n];
209        for (int i = n-1; i >= 0; i--) {
210            for (int j = n-i; j > 0; j--) {
211                coefficients[j] = coefficients[j-1] - c[i] * coefficients[j];
212            }
213            coefficients[0] = a[i] - c[i] * coefficients[0];
214        }
215
216        coefficientsComputed = true;
217    }
218
219    /**
220     * Verifies that the input arrays are valid.
221     * <p>
222     * The centers must be distinct for interpolation purposes, but not
223     * for general use. Thus it is not verified here.</p>
224     *
225     * @param a the coefficients in Newton form formula
226     * @param c the centers
227     * @throws NullArgumentException if any argument is {@code null}.
228     * @throws NoDataException if any array has zero length.
229     * @throws DimensionMismatchException if the size difference between
230     * {@code a} and {@code c} is not equal to 1.
231     * @see org.apache.commons.math4.analysis.interpolation.DividedDifferenceInterpolator#computeDividedDifference(double[],
232     * double[])
233     */
234    protected static void verifyInputArray(double a[], double c[])
235        throws NullArgumentException, NoDataException, DimensionMismatchException {
236        MathUtils.checkNotNull(a);
237        MathUtils.checkNotNull(c);
238        if (a.length == 0 || c.length == 0) {
239            throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY);
240        }
241        if (a.length != c.length + 1) {
242            throw new DimensionMismatchException(LocalizedFormats.ARRAY_SIZES_SHOULD_HAVE_DIFFERENCE_1,
243                                                 a.length, c.length);
244        }
245    }
246
247}