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