Classes in this File | Line Coverage | Branch Coverage | Complexity | ||||||||
PolynomialFunctionNewtonForm |
|
| 2.0;2 |
1 | /* |
|
2 | * Copyright 2003-2005 The Apache Software Foundation. |
|
3 | * |
|
4 | * Licensed under the Apache License, Version 2.0 (the "License"); |
|
5 | * you may not use this file except in compliance with the License. |
|
6 | * You may obtain a copy of the License at |
|
7 | * |
|
8 | * http://www.apache.org/licenses/LICENSE-2.0 |
|
9 | * |
|
10 | * Unless required by applicable law or agreed to in writing, software |
|
11 | * distributed under the License is distributed on an "AS IS" BASIS, |
|
12 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
13 | * See the License for the specific language governing permissions and |
|
14 | * limitations under the License. |
|
15 | */ |
|
16 | package org.apache.commons.math.analysis; |
|
17 | ||
18 | import java.io.Serializable; |
|
19 | import org.apache.commons.math.FunctionEvaluationException; |
|
20 | ||
21 | /** |
|
22 | * Implements the representation of a real polynomial function in |
|
23 | * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>, |
|
24 | * ISBN 0070124477, chapter 2. |
|
25 | * <p> |
|
26 | * The formula of polynomial in Newton form is |
|
27 | * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + |
|
28 | * a[n](x-c[0])(x-c[1])...(x-c[n-1]) |
|
29 | * Note that the length of a[] is one more than the length of c[] |
|
30 | * |
|
31 | * @version $Revision$ $Date$ |
|
32 | */ |
|
33 | public class PolynomialFunctionNewtonForm implements UnivariateRealFunction, |
|
34 | Serializable { |
|
35 | ||
36 | /** serializable version identifier */ |
|
37 | static final long serialVersionUID = -3353896576191389897L; |
|
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. |
|
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 | PolynomialFunctionNewtonForm(double a[], double c[]) throws |
|
70 | 14 | IllegalArgumentException { |
71 | ||
72 | 14 | verifyInputArray(a, c); |
73 | 10 | this.a = new double[a.length]; |
74 | 10 | this.c = new double[c.length]; |
75 | 10 | System.arraycopy(a, 0, this.a, 0, a.length); |
76 | 10 | System.arraycopy(c, 0, this.c, 0, c.length); |
77 | 10 | coefficientsComputed = false; |
78 | 10 | } |
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 | 28 | 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 | 12 | 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. |
|
105 | * |
|
106 | * @return a fresh copy of coefficients in Newton form formula |
|
107 | */ |
|
108 | public double[] getNewtonCoefficients() { |
|
109 | 0 | double[] out = new double[a.length]; |
110 | 0 | System.arraycopy(a, 0, out, 0, a.length); |
111 | 0 | 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. |
|
118 | * |
|
119 | * @return a fresh copy of the centers array |
|
120 | */ |
|
121 | public double[] getCenters() { |
|
122 | 0 | double[] out = new double[c.length]; |
123 | 0 | System.arraycopy(c, 0, out, 0, c.length); |
124 | 0 | 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. |
|
131 | * |
|
132 | * @return a fresh copy of the coefficients array |
|
133 | */ |
|
134 | public double[] getCoefficients() { |
|
135 | 6 | if (!coefficientsComputed) { |
136 | 6 | computeCoefficients(); |
137 | } |
|
138 | 6 | double[] out = new double[coefficients.length]; |
139 | 6 | System.arraycopy(coefficients, 0, out, 0, coefficients.length); |
140 | 6 | 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 | 28 | verifyInputArray(a, c); |
159 | ||
160 | 28 | int n = c.length; |
161 | 28 | double value = a[n]; |
162 | 120 | for (int i = n-1; i >= 0; i--) { |
163 | 92 | value = a[i] + (z - c[i]) * value; |
164 | } |
|
165 | ||
166 | 28 | 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 | 6 | int i, j, n = degree(); |
175 | ||
176 | 6 | coefficients = new double[n+1]; |
177 | 28 | for (i = 0; i <= n; i++) { |
178 | 22 | coefficients[i] = 0.0; |
179 | } |
|
180 | ||
181 | 6 | coefficients[0] = a[n]; |
182 | 22 | for (i = n-1; i >= 0; i--) { |
183 | 54 | for (j = n-i; j > 0; j--) { |
184 | 38 | coefficients[j] = coefficients[j-1] - c[i] * coefficients[j]; |
185 | } |
|
186 | 16 | coefficients[0] = a[i] - c[i] * coefficients[0]; |
187 | } |
|
188 | ||
189 | 6 | coefficientsComputed = true; |
190 | 6 | } |
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. |
|
197 | * |
|
198 | * @throws IllegalArgumentException if not valid |
|
199 | * @see DividedDifferenceInterpolator#computeDividedDifference(double[], |
|
200 | * double[]) |
|
201 | */ |
|
202 | protected static void verifyInputArray(double a[], double c[]) throws |
|
203 | IllegalArgumentException { |
|
204 | ||
205 | 42 | if (a.length < 1 || c.length < 1) { |
206 | 0 | throw new IllegalArgumentException |
207 | ("Input arrays must not be empty."); |
|
208 | } |
|
209 | 42 | if (a.length != c.length + 1) { |
210 | 4 | throw new IllegalArgumentException |
211 | ("Bad input array sizes, should have difference 1."); |
|
212 | } |
|
213 | 38 | } |
214 | } |