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