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