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    }