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