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    }