View Javadoc

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 }