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.MathRuntimeException; 020 import org.apache.commons.math.analysis.UnivariateRealFunction; 021 import org.apache.commons.math.FunctionEvaluationException; 022 import org.apache.commons.math.exception.util.LocalizedFormats; 023 024 /** 025 * Implements the representation of a real polynomial function in 026 * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>, 027 * ISBN 0070124477, chapter 2. 028 * <p> 029 * The formula of polynomial in Newton form is 030 * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + 031 * a[n](x-c[0])(x-c[1])...(x-c[n-1]) 032 * Note that the length of a[] is one more than the length of c[]</p> 033 * 034 * @version $Revision: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 f??vr. 2011) $ 035 * @since 1.2 036 */ 037 public class PolynomialFunctionNewtonForm implements UnivariateRealFunction { 038 039 /** 040 * The coefficients of the polynomial, ordered by degree -- i.e. 041 * coefficients[0] is the constant term and coefficients[n] is the 042 * coefficient of x^n where n is the degree of the polynomial. 043 */ 044 private double coefficients[]; 045 046 /** 047 * Centers of the Newton polynomial. 048 */ 049 private final double c[]; 050 051 /** 052 * When all c[i] = 0, a[] becomes normal polynomial coefficients, 053 * i.e. a[i] = coefficients[i]. 054 */ 055 private final double a[]; 056 057 /** 058 * Whether the polynomial coefficients are available. 059 */ 060 private boolean coefficientsComputed; 061 062 /** 063 * Construct a Newton polynomial with the given a[] and c[]. The order of 064 * centers are important in that if c[] shuffle, then values of a[] would 065 * completely change, not just a permutation of old a[]. 066 * <p> 067 * The constructor makes copy of the input arrays and assigns them.</p> 068 * 069 * @param a the coefficients in Newton form formula 070 * @param c the centers 071 * @throws IllegalArgumentException if input arrays are not valid 072 */ 073 public PolynomialFunctionNewtonForm(double a[], double c[]) 074 throws IllegalArgumentException { 075 076 verifyInputArray(a, c); 077 this.a = new double[a.length]; 078 this.c = new double[c.length]; 079 System.arraycopy(a, 0, this.a, 0, a.length); 080 System.arraycopy(c, 0, this.c, 0, c.length); 081 coefficientsComputed = false; 082 } 083 084 /** 085 * Calculate the function value at the given point. 086 * 087 * @param z the point at which the function value is to be computed 088 * @return the function value 089 * @throws FunctionEvaluationException if a runtime error occurs 090 * @see UnivariateRealFunction#value(double) 091 */ 092 public double value(double z) throws FunctionEvaluationException { 093 return evaluate(a, c, z); 094 } 095 096 /** 097 * Returns the degree of the polynomial. 098 * 099 * @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 }