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 */
17 package org.apache.commons.math.analysis.polynomials;
18
19 import org.apache.commons.math.FunctionEvaluationException;
20 import org.apache.commons.math.MathRuntimeException;
21 import org.apache.commons.math.analysis.UnivariateRealFunction;
22 import org.apache.commons.math.analysis.interpolation.DividedDifferenceInterpolator;
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: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
35 * @since 1.2
36 */
37 public 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 * Members of c[] are called centers of the Newton polynomial.
48 * When all c[i] = 0, a[] becomes normal polynomial coefficients,
49 * i.e. a[i] = coefficients[i].
50 */
51 private double a[], c[];
52
53 /**
54 * Whether the polynomial coefficients are available.
55 */
56 private boolean coefficientsComputed;
57
58 /**
59 * Construct a Newton polynomial with the given a[] and c[]. The order of
60 * centers are important in that if c[] shuffle, then values of a[] would
61 * completely change, not just a permutation of old a[].
62 * <p>
63 * The constructor makes copy of the input arrays and assigns them.</p>
64 *
65 * @param a the coefficients in Newton form formula
66 * @param c the centers
67 * @throws IllegalArgumentException if input arrays are not valid
68 */
69 public PolynomialFunctionNewtonForm(double a[], double c[])
70 throws IllegalArgumentException {
71
72 verifyInputArray(a, c);
73 this.a = new double[a.length];
74 this.c = new double[c.length];
75 System.arraycopy(a, 0, this.a, 0, a.length);
76 System.arraycopy(c, 0, this.c, 0, c.length);
77 coefficientsComputed = false;
78 }
79
80 /**
81 * Calculate the function value at the given point.
82 *
83 * @param z the point at which the function value is to be computed
84 * @return the function value
85 * @throws FunctionEvaluationException if a runtime error occurs
86 * @see UnivariateRealFunction#value(double)
87 */
88 public double value(double z) throws FunctionEvaluationException {
89 return evaluate(a, c, z);
90 }
91
92 /**
93 * Returns the degree of the polynomial.
94 *
95 * @return the degree of the polynomial
96 */
97 public int degree() {
98 return c.length;
99 }
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 }