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 }