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