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