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; 18 19 import java.io.Serializable; 20 import org.apache.commons.math.FunctionEvaluationException; 21 22 /** 23 * Implements the representation of a real polynomial function in 24 * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>, 25 * ISBN 0070124477, chapter 2. 26 * <p> 27 * The formula of polynomial in Newton form is 28 * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + 29 * a[n](x-c[0])(x-c[1])...(x-c[n-1]) 30 * Note that the length of a[] is one more than the length of c[]</p> 31 * 32 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $ 33 * @since 1.2 34 */ 35 public class PolynomialFunctionNewtonForm implements UnivariateRealFunction, 36 Serializable { 37 38 /** serializable version identifier */ 39 static final long serialVersionUID = -3353896576191389897L; 40 41 /** 42 * The coefficients of the polynomial, ordered by degree -- i.e. 43 * coefficients[0] is the constant term and coefficients[n] is the 44 * coefficient of x^n where n is the degree of the polynomial. 45 */ 46 private double coefficients[]; 47 48 /** 49 * Members of c[] are called centers of the Newton polynomial. 50 * When all c[i] = 0, a[] becomes normal polynomial coefficients, 51 * i.e. a[i] = coefficients[i]. 52 */ 53 private double a[], c[]; 54 55 /** 56 * Whether the polynomial coefficients are available. 57 */ 58 private boolean coefficientsComputed; 59 60 /** 61 * Construct a Newton polynomial with the given a[] and c[]. The order of 62 * centers are important in that if c[] shuffle, then values of a[] would 63 * completely change, not just a permutation of old a[]. 64 * <p> 65 * The constructor makes copy of the input arrays and assigns them.</p> 66 * 67 * @param a the coefficients in Newton form formula 68 * @param c the centers 69 * @throws IllegalArgumentException if input arrays are not valid 70 */ 71 PolynomialFunctionNewtonForm(double a[], double c[]) throws 72 IllegalArgumentException { 73 74 verifyInputArray(a, c); 75 this.a = new double[a.length]; 76 this.c = new double[c.length]; 77 System.arraycopy(a, 0, this.a, 0, a.length); 78 System.arraycopy(c, 0, this.c, 0, c.length); 79 coefficientsComputed = false; 80 } 81 82 /** 83 * Calculate the function value at the given point. 84 * 85 * @param z the point at which the function value is to be computed 86 * @return the function value 87 * @throws FunctionEvaluationException if a runtime error occurs 88 * @see UnivariateRealFunction#value(double) 89 */ 90 public double value(double z) throws FunctionEvaluationException { 91 return evaluate(a, c, z); 92 } 93 94 /** 95 * Returns the degree of the polynomial. 96 * 97 * @return the degree of the polynomial 98 */ 99 public int degree() { 100 return c.length; 101 } 102 103 /** 104 * Returns a copy of coefficients in Newton form formula. 105 * <p> 106 * Changes made to the returned copy will not affect the polynomial.</p> 107 * 108 * @return a fresh copy of coefficients in Newton form formula 109 */ 110 public double[] getNewtonCoefficients() { 111 double[] out = new double[a.length]; 112 System.arraycopy(a, 0, out, 0, a.length); 113 return out; 114 } 115 116 /** 117 * Returns a copy of the centers array. 118 * <p> 119 * Changes made to the returned copy will not affect the polynomial.</p> 120 * 121 * @return a fresh copy of the centers array 122 */ 123 public double[] getCenters() { 124 double[] out = new double[c.length]; 125 System.arraycopy(c, 0, out, 0, c.length); 126 return out; 127 } 128 129 /** 130 * Returns a copy of the coefficients array. 131 * <p> 132 * Changes made to the returned copy will not affect the polynomial.</p> 133 * 134 * @return a fresh copy of the coefficients array 135 */ 136 public double[] getCoefficients() { 137 if (!coefficientsComputed) { 138 computeCoefficients(); 139 } 140 double[] out = new double[coefficients.length]; 141 System.arraycopy(coefficients, 0, out, 0, coefficients.length); 142 return out; 143 } 144 145 /** 146 * Evaluate the Newton polynomial using nested multiplication. It is 147 * also called <a href="http://mathworld.wolfram.com/HornersRule.html"> 148 * Horner's Rule</a> and takes O(N) time. 149 * 150 * @param a the coefficients in Newton form formula 151 * @param c the centers 152 * @param z the point at which the function value is to be computed 153 * @return the function value 154 * @throws FunctionEvaluationException if a runtime error occurs 155 * @throws IllegalArgumentException if inputs are not valid 156 */ 157 public static double evaluate(double a[], double c[], double z) throws 158 FunctionEvaluationException, IllegalArgumentException { 159 160 verifyInputArray(a, c); 161 162 int n = c.length; 163 double value = a[n]; 164 for (int i = n-1; i >= 0; i--) { 165 value = a[i] + (z - c[i]) * value; 166 } 167 168 return value; 169 } 170 171 /** 172 * Calculate the normal polynomial coefficients given the Newton form. 173 * It also uses nested multiplication but takes O(N^2) time. 174 */ 175 protected void computeCoefficients() { 176 int i, j, n = degree(); 177 178 coefficients = new double[n+1]; 179 for (i = 0; i <= n; i++) { 180 coefficients[i] = 0.0; 181 } 182 183 coefficients[0] = a[n]; 184 for (i = n-1; i >= 0; i--) { 185 for (j = n-i; j > 0; j--) { 186 coefficients[j] = coefficients[j-1] - c[i] * coefficients[j]; 187 } 188 coefficients[0] = a[i] - c[i] * coefficients[0]; 189 } 190 191 coefficientsComputed = true; 192 } 193 194 /** 195 * Verifies that the input arrays are valid. 196 * <p> 197 * The centers must be distinct for interpolation purposes, but not 198 * for general use. Thus it is not verified here.</p> 199 * 200 * @param a the coefficients in Newton form formula 201 * @param c the centers 202 * @throws IllegalArgumentException if not valid 203 * @see DividedDifferenceInterpolator#computeDividedDifference(double[], 204 * double[]) 205 */ 206 protected static void verifyInputArray(double a[], double c[]) throws 207 IllegalArgumentException { 208 209 if (a.length < 1 || c.length < 1) { 210 throw new IllegalArgumentException 211 ("Input arrays must not be empty."); 212 } 213 if (a.length != c.length + 1) { 214 throw new IllegalArgumentException 215 ("Bad input array sizes, should have difference 1."); 216 } 217 } 218 }