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