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.FunctionEvaluationException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.analysis.UnivariateRealFunction; 022 023 /** 024 * Implements the representation of a real polynomial function in 025 * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>, 026 * ISBN 0070124477, chapter 2. 027 * <p> 028 * The formula of polynomial in Newton form is 029 * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + 030 * a[n](x-c[0])(x-c[1])...(x-c[n-1]) 031 * Note that the length of a[] is one more than the length of c[]</p> 032 * 033 * @version $Revision: 922708 $ $Date: 2010-03-13 20:15:47 -0500 (Sat, 13 Mar 2010) $ 034 * @since 1.2 035 */ 036 public class PolynomialFunctionNewtonForm implements UnivariateRealFunction { 037 038 /** 039 * The coefficients of the polynomial, ordered by degree -- i.e. 040 * coefficients[0] is the constant term and coefficients[n] is the 041 * coefficient of x^n where n is the degree of the polynomial. 042 */ 043 private double coefficients[]; 044 045 /** 046 * Centers of the Newton polynomial. 047 */ 048 private final double c[]; 049 050 /** 051 * When all c[i] = 0, a[] becomes normal polynomial coefficients, 052 * i.e. a[i] = coefficients[i]. 053 */ 054 private final double a[]; 055 056 /** 057 * Whether the polynomial coefficients are available. 058 */ 059 private boolean coefficientsComputed; 060 061 /** 062 * Construct a Newton polynomial with the given a[] and c[]. The order of 063 * centers are important in that if c[] shuffle, then values of a[] would 064 * completely change, not just a permutation of old a[]. 065 * <p> 066 * The constructor makes copy of the input arrays and assigns them.</p> 067 * 068 * @param a the coefficients in Newton form formula 069 * @param c the centers 070 * @throws IllegalArgumentException if input arrays are not valid 071 */ 072 public PolynomialFunctionNewtonForm(double a[], double c[]) 073 throws IllegalArgumentException { 074 075 verifyInputArray(a, c); 076 this.a = new double[a.length]; 077 this.c = new double[c.length]; 078 System.arraycopy(a, 0, this.a, 0, a.length); 079 System.arraycopy(c, 0, this.c, 0, c.length); 080 coefficientsComputed = false; 081 } 082 083 /** 084 * Calculate the function value at the given point. 085 * 086 * @param z the point at which the function value is to be computed 087 * @return the function value 088 * @throws FunctionEvaluationException if a runtime error occurs 089 * @see UnivariateRealFunction#value(double) 090 */ 091 public double value(double z) throws FunctionEvaluationException { 092 return evaluate(a, c, z); 093 } 094 095 /** 096 * Returns the degree of the polynomial. 097 * 098 * @return the degree of the polynomial 099 */ 100 public int degree() { 101 return c.length; 102 } 103 104 /** 105 * Returns a copy of coefficients in Newton form formula. 106 * <p> 107 * Changes made to the returned copy will not affect the polynomial.</p> 108 * 109 * @return a fresh copy of coefficients in Newton form formula 110 */ 111 public double[] getNewtonCoefficients() { 112 double[] out = new double[a.length]; 113 System.arraycopy(a, 0, out, 0, a.length); 114 return out; 115 } 116 117 /** 118 * Returns a copy of the centers array. 119 * <p> 120 * Changes made to the returned copy will not affect the polynomial.</p> 121 * 122 * @return a fresh copy of the centers array 123 */ 124 public double[] getCenters() { 125 double[] out = new double[c.length]; 126 System.arraycopy(c, 0, out, 0, c.length); 127 return out; 128 } 129 130 /** 131 * Returns a copy of the coefficients array. 132 * <p> 133 * Changes made to the returned copy will not affect the polynomial.</p> 134 * 135 * @return a fresh copy of the coefficients array 136 */ 137 public double[] getCoefficients() { 138 if (!coefficientsComputed) { 139 computeCoefficients(); 140 } 141 double[] out = new double[coefficients.length]; 142 System.arraycopy(coefficients, 0, out, 0, coefficients.length); 143 return out; 144 } 145 146 /** 147 * Evaluate the Newton polynomial using nested multiplication. It is 148 * also called <a href="http://mathworld.wolfram.com/HornersRule.html"> 149 * Horner's Rule</a> and takes O(N) time. 150 * 151 * @param a the coefficients in Newton form formula 152 * @param c the centers 153 * @param z the point at which the function value is to be computed 154 * @return the function value 155 * @throws FunctionEvaluationException if a runtime error occurs 156 * @throws IllegalArgumentException if inputs are not valid 157 */ 158 public static double evaluate(double a[], double c[], double z) throws 159 FunctionEvaluationException, IllegalArgumentException { 160 161 verifyInputArray(a, c); 162 163 int n = c.length; 164 double value = a[n]; 165 for (int i = n-1; i >= 0; i--) { 166 value = a[i] + (z - c[i]) * value; 167 } 168 169 return value; 170 } 171 172 /** 173 * Calculate the normal polynomial coefficients given the Newton form. 174 * It also uses nested multiplication but takes O(N^2) time. 175 */ 176 protected void computeCoefficients() { 177 final int n = degree(); 178 179 coefficients = new double[n+1]; 180 for (int i = 0; i <= n; i++) { 181 coefficients[i] = 0.0; 182 } 183 184 coefficients[0] = a[n]; 185 for (int i = n-1; i >= 0; i--) { 186 for (int j = n-i; j > 0; j--) { 187 coefficients[j] = coefficients[j-1] - c[i] * coefficients[j]; 188 } 189 coefficients[0] = a[i] - c[i] * coefficients[0]; 190 } 191 192 coefficientsComputed = true; 193 } 194 195 /** 196 * Verifies that the input arrays are valid. 197 * <p> 198 * The centers must be distinct for interpolation purposes, but not 199 * for general use. Thus it is not verified here.</p> 200 * 201 * @param a the coefficients in Newton form formula 202 * @param c the centers 203 * @throws IllegalArgumentException if not valid 204 * @see org.apache.commons.math.analysis.interpolation.DividedDifferenceInterpolator#computeDividedDifference(double[], 205 * double[]) 206 */ 207 protected static void verifyInputArray(double a[], double c[]) throws 208 IllegalArgumentException { 209 210 if (a.length < 1 || c.length < 1) { 211 throw MathRuntimeException.createIllegalArgumentException( 212 "empty polynomials coefficients array"); 213 } 214 if (a.length != c.length + 1) { 215 throw MathRuntimeException.createIllegalArgumentException( 216 "array sizes should have difference 1 ({0} != {1} + 1)", 217 a.length, c.length); 218 } 219 } 220 }