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.MathRuntimeException; 20 import org.apache.commons.math.analysis.UnivariateRealFunction; 21 import org.apache.commons.math.FunctionEvaluationException; 22 import org.apache.commons.math.exception.util.LocalizedFormats; 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: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 fvr. 2011) $ 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 * Centers of the Newton polynomial. 48 */ 49 private final double c[]; 50 51 /** 52 * When all c[i] = 0, a[] becomes normal polynomial coefficients, 53 * i.e. a[i] = coefficients[i]. 54 */ 55 private final double a[]; 56 57 /** 58 * Whether the polynomial coefficients are available. 59 */ 60 private boolean coefficientsComputed; 61 62 /** 63 * Construct a Newton polynomial with the given a[] and c[]. The order of 64 * centers are important in that if c[] shuffle, then values of a[] would 65 * completely change, not just a permutation of old a[]. 66 * <p> 67 * The constructor makes copy of the input arrays and assigns them.</p> 68 * 69 * @param a the coefficients in Newton form formula 70 * @param c the centers 71 * @throws IllegalArgumentException if input arrays are not valid 72 */ 73 public PolynomialFunctionNewtonForm(double a[], double c[]) 74 throws IllegalArgumentException { 75 76 verifyInputArray(a, c); 77 this.a = new double[a.length]; 78 this.c = new double[c.length]; 79 System.arraycopy(a, 0, this.a, 0, a.length); 80 System.arraycopy(c, 0, this.c, 0, c.length); 81 coefficientsComputed = false; 82 } 83 84 /** 85 * Calculate the function value at the given point. 86 * 87 * @param z the point at which the function value is to be computed 88 * @return the function value 89 * @throws FunctionEvaluationException if a runtime error occurs 90 * @see UnivariateRealFunction#value(double) 91 */ 92 public double value(double z) throws FunctionEvaluationException { 93 return evaluate(a, c, z); 94 } 95 96 /** 97 * Returns the degree of the polynomial. 98 * 99 * @return the degree of the polynomial 100 */ 101 public int degree() { 102 return c.length; 103 } 104 105 /** 106 * Returns a copy of coefficients in Newton form formula. 107 * <p> 108 * Changes made to the returned copy will not affect the polynomial.</p> 109 * 110 * @return a fresh copy of coefficients in Newton form formula 111 */ 112 public double[] getNewtonCoefficients() { 113 double[] out = new double[a.length]; 114 System.arraycopy(a, 0, out, 0, a.length); 115 return out; 116 } 117 118 /** 119 * Returns a copy of the centers array. 120 * <p> 121 * Changes made to the returned copy will not affect the polynomial.</p> 122 * 123 * @return a fresh copy of the centers array 124 */ 125 public double[] getCenters() { 126 double[] out = new double[c.length]; 127 System.arraycopy(c, 0, out, 0, c.length); 128 return out; 129 } 130 131 /** 132 * Returns a copy of the coefficients array. 133 * <p> 134 * Changes made to the returned copy will not affect the polynomial.</p> 135 * 136 * @return a fresh copy of the coefficients array 137 */ 138 public double[] getCoefficients() { 139 if (!coefficientsComputed) { 140 computeCoefficients(); 141 } 142 double[] out = new double[coefficients.length]; 143 System.arraycopy(coefficients, 0, out, 0, coefficients.length); 144 return out; 145 } 146 147 /** 148 * Evaluate the Newton polynomial using nested multiplication. It is 149 * also called <a href="http://mathworld.wolfram.com/HornersRule.html"> 150 * Horner's Rule</a> and takes O(N) time. 151 * 152 * @param a the coefficients in Newton form formula 153 * @param c the centers 154 * @param z the point at which the function value is to be computed 155 * @return the function value 156 * @throws FunctionEvaluationException if a runtime error occurs 157 * @throws IllegalArgumentException if inputs are not valid 158 */ 159 public static double evaluate(double a[], double c[], double z) 160 throws FunctionEvaluationException, IllegalArgumentException { 161 162 verifyInputArray(a, c); 163 164 int n = c.length; 165 double value = a[n]; 166 for (int i = n-1; i >= 0; i--) { 167 value = a[i] + (z - c[i]) * value; 168 } 169 170 return value; 171 } 172 173 /** 174 * Calculate the normal polynomial coefficients given the Newton form. 175 * It also uses nested multiplication but takes O(N^2) time. 176 */ 177 protected void computeCoefficients() { 178 final int n = degree(); 179 180 coefficients = new double[n+1]; 181 for (int i = 0; i <= n; i++) { 182 coefficients[i] = 0.0; 183 } 184 185 coefficients[0] = a[n]; 186 for (int i = n-1; i >= 0; i--) { 187 for (int j = n-i; j > 0; j--) { 188 coefficients[j] = coefficients[j-1] - c[i] * coefficients[j]; 189 } 190 coefficients[0] = a[i] - c[i] * coefficients[0]; 191 } 192 193 coefficientsComputed = true; 194 } 195 196 /** 197 * Verifies that the input arrays are valid. 198 * <p> 199 * The centers must be distinct for interpolation purposes, but not 200 * for general use. Thus it is not verified here.</p> 201 * 202 * @param a the coefficients in Newton form formula 203 * @param c the centers 204 * @throws IllegalArgumentException if not valid 205 * @see org.apache.commons.math.analysis.interpolation.DividedDifferenceInterpolator#computeDividedDifference(double[], 206 * double[]) 207 */ 208 protected static void verifyInputArray(double a[], double c[]) throws 209 IllegalArgumentException { 210 211 if (a.length < 1 || c.length < 1) { 212 throw MathRuntimeException.createIllegalArgumentException( 213 LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY); 214 } 215 if (a.length != c.length + 1) { 216 throw MathRuntimeException.createIllegalArgumentException( 217 LocalizedFormats.ARRAY_SIZES_SHOULD_HAVE_DIFFERENCE_1, 218 a.length, c.length); 219 } 220 } 221 } 222