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 import org.apache.commons.math.analysis.interpolation.DividedDifferenceInterpolator; 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: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 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 * Members of c[] are called centers of the Newton polynomial. 048 * When all c[i] = 0, a[] becomes normal polynomial coefficients, 049 * i.e. a[i] = coefficients[i]. 050 */ 051 private double a[], c[]; 052 053 /** 054 * Whether the polynomial coefficients are available. 055 */ 056 private boolean coefficientsComputed; 057 058 /** 059 * Construct a Newton polynomial with the given a[] and c[]. The order of 060 * centers are important in that if c[] shuffle, then values of a[] would 061 * completely change, not just a permutation of old a[]. 062 * <p> 063 * The constructor makes copy of the input arrays and assigns them.</p> 064 * 065 * @param a the coefficients in Newton form formula 066 * @param c the centers 067 * @throws IllegalArgumentException if input arrays are not valid 068 */ 069 public PolynomialFunctionNewtonForm(double a[], double c[]) 070 throws IllegalArgumentException { 071 072 verifyInputArray(a, c); 073 this.a = new double[a.length]; 074 this.c = new double[c.length]; 075 System.arraycopy(a, 0, this.a, 0, a.length); 076 System.arraycopy(c, 0, this.c, 0, c.length); 077 coefficientsComputed = false; 078 } 079 080 /** 081 * Calculate the function value at the given point. 082 * 083 * @param z the point at which the function value is to be computed 084 * @return the function value 085 * @throws FunctionEvaluationException if a runtime error occurs 086 * @see UnivariateRealFunction#value(double) 087 */ 088 public double value(double z) throws FunctionEvaluationException { 089 return evaluate(a, c, z); 090 } 091 092 /** 093 * Returns the degree of the polynomial. 094 * 095 * @return the degree of the polynomial 096 */ 097 public int degree() { 098 return c.length; 099 } 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 }