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