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