1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.analysis.polynomials;
18
19import org.apache.commons.math.MathRuntimeException;
20import org.apache.commons.math.analysis.UnivariateRealFunction;
21import org.apache.commons.math.FunctionEvaluationException;
22import org.apache.commons.math.exception.util.LocalizedFormats;
23
24/**
25 * Implements the representation of a real polynomial function in
26 * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>,
27 * ISBN 0070124477, chapter 2.
28 * <p>
29 * The formula of polynomial in Newton form is
30 *     p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
31 *            a[n](x-c[0])(x-c[1])...(x-c[n-1])
32 * Note that the length of a[] is one more than the length of c[]</p>
33 *
34 * @version $Revision: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 févr. 2011) $
35 * @since 1.2
36 */
37public class PolynomialFunctionNewtonForm implements UnivariateRealFunction {
38
39    /**
40     * The coefficients of the polynomial, ordered by degree -- i.e.
41     * coefficients[0] is the constant term and coefficients[n] is the
42     * coefficient of x^n where n is the degree of the polynomial.
43     */
44    private double coefficients[];
45
46    /**
47     * Centers of the Newton polynomial.
48     */
49    private final double c[];
50
51    /**
52     * When all c[i] = 0, a[] becomes normal polynomial coefficients,
53     * i.e. a[i] = coefficients[i].
54     */
55    private final double a[];
56
57    /**
58     * Whether the polynomial coefficients are available.
59     */
60    private boolean coefficientsComputed;
61
62    /**
63     * Construct a Newton polynomial with the given a[] and c[]. The order of
64     * centers are important in that if c[] shuffle, then values of a[] would
65     * completely change, not just a permutation of old a[].
66     * <p>
67     * The constructor makes copy of the input arrays and assigns them.</p>
68     *
69     * @param a the coefficients in Newton form formula
70     * @param c the centers
71     * @throws IllegalArgumentException if input arrays are not valid
72     */
73    public PolynomialFunctionNewtonForm(double a[], double c[])
74        throws IllegalArgumentException {
75
76        verifyInputArray(a, c);
77        this.a = new double[a.length];
78        this.c = new double[c.length];
79        System.arraycopy(a, 0, this.a, 0, a.length);
80        System.arraycopy(c, 0, this.c, 0, c.length);
81        coefficientsComputed = false;
82    }
83
84    /**
85     * Calculate the function value at the given point.
86     *
87     * @param z the point at which the function value is to be computed
88     * @return the function value
89     * @throws FunctionEvaluationException if a runtime error occurs
90     * @see UnivariateRealFunction#value(double)
91     */
92    public double value(double z) throws FunctionEvaluationException {
93       return evaluate(a, c, z);
94    }
95
96    /**
97     * Returns the degree of the polynomial.
98     *
99     * @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}
222