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 java.util.ArrayList; 020 import java.util.HashMap; 021 import java.util.List; 022 import java.util.Map; 023 024 import org.apache.commons.math3.fraction.BigFraction; 025 import org.apache.commons.math3.util.ArithmeticUtils; 026 import org.apache.commons.math3.util.FastMath; 027 028 /** 029 * A collection of static methods that operate on or return polynomials. 030 * 031 * @version $Id: PolynomialsUtils.java 1364387 2012-07-22 18:14:11Z tn $ 032 * @since 2.0 033 */ 034 public class PolynomialsUtils { 035 036 /** Coefficients for Chebyshev polynomials. */ 037 private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS; 038 039 /** Coefficients for Hermite polynomials. */ 040 private static final List<BigFraction> HERMITE_COEFFICIENTS; 041 042 /** Coefficients for Laguerre polynomials. */ 043 private static final List<BigFraction> LAGUERRE_COEFFICIENTS; 044 045 /** Coefficients for Legendre polynomials. */ 046 private static final List<BigFraction> LEGENDRE_COEFFICIENTS; 047 048 /** Coefficients for Jacobi polynomials. */ 049 private static final Map<JacobiKey, List<BigFraction>> JACOBI_COEFFICIENTS; 050 051 static { 052 053 // initialize recurrence for Chebyshev polynomials 054 // T0(X) = 1, T1(X) = 0 + 1 * X 055 CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>(); 056 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE); 057 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO); 058 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE); 059 060 // initialize recurrence for Hermite polynomials 061 // H0(X) = 1, H1(X) = 0 + 2 * X 062 HERMITE_COEFFICIENTS = new ArrayList<BigFraction>(); 063 HERMITE_COEFFICIENTS.add(BigFraction.ONE); 064 HERMITE_COEFFICIENTS.add(BigFraction.ZERO); 065 HERMITE_COEFFICIENTS.add(BigFraction.TWO); 066 067 // initialize recurrence for Laguerre polynomials 068 // L0(X) = 1, L1(X) = 1 - 1 * X 069 LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>(); 070 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE); 071 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE); 072 LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE); 073 074 // initialize recurrence for Legendre polynomials 075 // P0(X) = 1, P1(X) = 0 + 1 * X 076 LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>(); 077 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE); 078 LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO); 079 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE); 080 081 // initialize map for Jacobi polynomials 082 JACOBI_COEFFICIENTS = new HashMap<JacobiKey, List<BigFraction>>(); 083 084 } 085 086 /** 087 * Private constructor, to prevent instantiation. 088 */ 089 private PolynomialsUtils() { 090 } 091 092 /** 093 * Create a Chebyshev polynomial of the first kind. 094 * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev 095 * polynomials of the first kind</a> are orthogonal polynomials. 096 * They can be defined by the following recurrence relations: 097 * <pre> 098 * T<sub>0</sub>(X) = 1 099 * T<sub>1</sub>(X) = X 100 * T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X) 101 * </pre></p> 102 * @param degree degree of the polynomial 103 * @return Chebyshev polynomial of specified degree 104 */ 105 public static PolynomialFunction createChebyshevPolynomial(final int degree) { 106 return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS, 107 new RecurrenceCoefficientsGenerator() { 108 private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE }; 109 /** {@inheritDoc} */ 110 public BigFraction[] generate(int k) { 111 return coeffs; 112 } 113 }); 114 } 115 116 /** 117 * Create a Hermite polynomial. 118 * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite 119 * polynomials</a> are orthogonal polynomials. 120 * They can be defined by the following recurrence relations: 121 * <pre> 122 * H<sub>0</sub>(X) = 1 123 * H<sub>1</sub>(X) = 2X 124 * H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X) 125 * </pre></p> 126 127 * @param degree degree of the polynomial 128 * @return Hermite polynomial of specified degree 129 */ 130 public static PolynomialFunction createHermitePolynomial(final int degree) { 131 return buildPolynomial(degree, HERMITE_COEFFICIENTS, 132 new RecurrenceCoefficientsGenerator() { 133 /** {@inheritDoc} */ 134 public BigFraction[] generate(int k) { 135 return new BigFraction[] { 136 BigFraction.ZERO, 137 BigFraction.TWO, 138 new BigFraction(2 * k)}; 139 } 140 }); 141 } 142 143 /** 144 * Create a Laguerre polynomial. 145 * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre 146 * polynomials</a> are orthogonal polynomials. 147 * They can be defined by the following recurrence relations: 148 * <pre> 149 * L<sub>0</sub>(X) = 1 150 * L<sub>1</sub>(X) = 1 - X 151 * (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X) 152 * </pre></p> 153 * @param degree degree of the polynomial 154 * @return Laguerre polynomial of specified degree 155 */ 156 public static PolynomialFunction createLaguerrePolynomial(final int degree) { 157 return buildPolynomial(degree, LAGUERRE_COEFFICIENTS, 158 new RecurrenceCoefficientsGenerator() { 159 /** {@inheritDoc} */ 160 public BigFraction[] generate(int k) { 161 final int kP1 = k + 1; 162 return new BigFraction[] { 163 new BigFraction(2 * k + 1, kP1), 164 new BigFraction(-1, kP1), 165 new BigFraction(k, kP1)}; 166 } 167 }); 168 } 169 170 /** 171 * Create a Legendre polynomial. 172 * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre 173 * polynomials</a> are orthogonal polynomials. 174 * They can be defined by the following recurrence relations: 175 * <pre> 176 * P<sub>0</sub>(X) = 1 177 * P<sub>1</sub>(X) = X 178 * (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X) 179 * </pre></p> 180 * @param degree degree of the polynomial 181 * @return Legendre polynomial of specified degree 182 */ 183 public static PolynomialFunction createLegendrePolynomial(final int degree) { 184 return buildPolynomial(degree, LEGENDRE_COEFFICIENTS, 185 new RecurrenceCoefficientsGenerator() { 186 /** {@inheritDoc} */ 187 public BigFraction[] generate(int k) { 188 final int kP1 = k + 1; 189 return new BigFraction[] { 190 BigFraction.ZERO, 191 new BigFraction(k + kP1, kP1), 192 new BigFraction(k, kP1)}; 193 } 194 }); 195 } 196 197 /** 198 * Create a Jacobi polynomial. 199 * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi 200 * polynomials</a> are orthogonal polynomials. 201 * They can be defined by the following recurrence relations: 202 * <pre> 203 * P<sub>0</sub><sup>vw</sup>(X) = 1 204 * P<sub>-1</sub><sup>vw</sup>(X) = 0 205 * 2k(k + v + w)(2k + v + w - 2) P<sub>k</sub><sup>vw</sup>(X) = 206 * (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) X + v<sup>2</sup> - w<sup>2</sup>] P<sub>k-1</sub><sup>vw</sup>(X) 207 * - 2(k + v - 1)(k + w - 1)(2k + v + w) P<sub>k-2</sub><sup>vw</sup>(X) 208 * </pre></p> 209 * @param degree degree of the polynomial 210 * @param v first exponent 211 * @param w second exponent 212 * @return Jacobi polynomial of specified degree 213 */ 214 public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) { 215 216 // select the appropriate list 217 final JacobiKey key = new JacobiKey(v, w); 218 219 if (!JACOBI_COEFFICIENTS.containsKey(key)) { 220 221 // allocate a new list for v, w 222 final List<BigFraction> list = new ArrayList<BigFraction>(); 223 JACOBI_COEFFICIENTS.put(key, list); 224 225 // Pv,w,0(x) = 1; 226 list.add(BigFraction.ONE); 227 228 // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2 229 list.add(new BigFraction(v - w, 2)); 230 list.add(new BigFraction(2 + v + w, 2)); 231 232 } 233 234 return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key), 235 new RecurrenceCoefficientsGenerator() { 236 /** {@inheritDoc} */ 237 public BigFraction[] generate(int k) { 238 k++; 239 final int kvw = k + v + w; 240 final int twoKvw = kvw + k; 241 final int twoKvwM1 = twoKvw - 1; 242 final int twoKvwM2 = twoKvw - 2; 243 final int den = 2 * k * kvw * twoKvwM2; 244 245 return new BigFraction[] { 246 new BigFraction(twoKvwM1 * (v * v - w * w), den), 247 new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den), 248 new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den) 249 }; 250 } 251 }); 252 253 } 254 255 /** Inner class for Jacobi polynomials keys. */ 256 private static class JacobiKey { 257 258 /** First exponent. */ 259 private final int v; 260 261 /** Second exponent. */ 262 private final int w; 263 264 /** Simple constructor. 265 * @param v first exponent 266 * @param w second exponent 267 */ 268 public JacobiKey(final int v, final int w) { 269 this.v = v; 270 this.w = w; 271 } 272 273 /** Get hash code. 274 * @return hash code 275 */ 276 @Override 277 public int hashCode() { 278 return (v << 16) ^ w; 279 } 280 281 /** Check if the instance represent the same key as another instance. 282 * @param key other key 283 * @return true if the instance and the other key refer to the same polynomial 284 */ 285 @Override 286 public boolean equals(final Object key) { 287 288 if ((key == null) || !(key instanceof JacobiKey)) { 289 return false; 290 } 291 292 final JacobiKey otherK = (JacobiKey) key; 293 return (v == otherK.v) && (w == otherK.w); 294 295 } 296 } 297 298 /** 299 * Compute the coefficients of the polynomial <code>P<sub>s</sub>(x)</code> 300 * whose values at point {@code x} will be the same as the those from the 301 * original polynomial <code>P(x)</code> when computed at {@code x + shift}. 302 * Thus, if <code>P(x) = Σ<sub>i</sub> a<sub>i</sub> x<sup>i</sup></code>, 303 * then 304 * <pre> 305 * <table> 306 * <tr> 307 * <td><code>P<sub>s</sub>(x)</td> 308 * <td>= Σ<sub>i</sub> b<sub>i</sub> x<sup>i</sup></code></td> 309 * </tr> 310 * <tr> 311 * <td></td> 312 * <td>= Σ<sub>i</sub> a<sub>i</sub> (x + shift)<sup>i</sup></code></td> 313 * </tr> 314 * </table> 315 * </pre> 316 * 317 * @param coefficients Coefficients of the original polynomial. 318 * @param shift Shift value. 319 * @return the coefficients <code>b<sub>i</sub></code> of the shifted 320 * polynomial. 321 */ 322 public static double[] shift(final double[] coefficients, 323 final double shift) { 324 final int dp1 = coefficients.length; 325 final double[] newCoefficients = new double[dp1]; 326 327 // Pascal triangle. 328 final int[][] coeff = new int[dp1][dp1]; 329 for (int i = 0; i < dp1; i++){ 330 for(int j = 0; j <= i; j++){ 331 coeff[i][j] = (int) ArithmeticUtils.binomialCoefficient(i, j); 332 } 333 } 334 335 // First polynomial coefficient. 336 for (int i = 0; i < dp1; i++){ 337 newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i); 338 } 339 340 // Superior order. 341 final int d = dp1 - 1; 342 for (int i = 0; i < d; i++) { 343 for (int j = i; j < d; j++){ 344 newCoefficients[i + 1] += coeff[j + 1][j - i] * 345 coefficients[j + 1] * FastMath.pow(shift, j - i); 346 } 347 } 348 349 return newCoefficients; 350 } 351 352 353 /** Get the coefficients array for a given degree. 354 * @param degree degree of the polynomial 355 * @param coefficients list where the computed coefficients are stored 356 * @param generator recurrence coefficients generator 357 * @return coefficients array 358 */ 359 private static PolynomialFunction buildPolynomial(final int degree, 360 final List<BigFraction> coefficients, 361 final RecurrenceCoefficientsGenerator generator) { 362 363 final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1; 364 synchronized (PolynomialsUtils.class) { 365 if (degree > maxDegree) { 366 computeUpToDegree(degree, maxDegree, generator, coefficients); 367 } 368 } 369 370 // coefficient for polynomial 0 is l [0] 371 // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1) 372 // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2) 373 // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3) 374 // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4) 375 // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5) 376 // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6) 377 // ... 378 final int start = degree * (degree + 1) / 2; 379 380 final double[] a = new double[degree + 1]; 381 for (int i = 0; i <= degree; ++i) { 382 a[i] = coefficients.get(start + i).doubleValue(); 383 } 384 385 // build the polynomial 386 return new PolynomialFunction(a); 387 388 } 389 390 /** Compute polynomial coefficients up to a given degree. 391 * @param degree maximal degree 392 * @param maxDegree current maximal degree 393 * @param generator recurrence coefficients generator 394 * @param coefficients list where the computed coefficients should be appended 395 */ 396 private static void computeUpToDegree(final int degree, final int maxDegree, 397 final RecurrenceCoefficientsGenerator generator, 398 final List<BigFraction> coefficients) { 399 400 int startK = (maxDegree - 1) * maxDegree / 2; 401 for (int k = maxDegree; k < degree; ++k) { 402 403 // start indices of two previous polynomials Pk(X) and Pk-1(X) 404 int startKm1 = startK; 405 startK += k; 406 407 // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X) 408 BigFraction[] ai = generator.generate(k); 409 410 BigFraction ck = coefficients.get(startK); 411 BigFraction ckm1 = coefficients.get(startKm1); 412 413 // degree 0 coefficient 414 coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2]))); 415 416 // degree 1 to degree k-1 coefficients 417 for (int i = 1; i < k; ++i) { 418 final BigFraction ckPrev = ck; 419 ck = coefficients.get(startK + i); 420 ckm1 = coefficients.get(startKm1 + i); 421 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2]))); 422 } 423 424 // degree k coefficient 425 final BigFraction ckPrev = ck; 426 ck = coefficients.get(startK + k); 427 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1]))); 428 429 // degree k+1 coefficient 430 coefficients.add(ck.multiply(ai[1])); 431 432 } 433 434 } 435 436 /** Interface for recurrence coefficients generation. */ 437 private interface RecurrenceCoefficientsGenerator { 438 /** 439 * Generate recurrence coefficients. 440 * @param k highest degree of the polynomials used in the recurrence 441 * @return an array of three coefficients such that 442 * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X) 443 */ 444 BigFraction[] generate(int k); 445 } 446 447 }