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="https://en.wikipedia.org/wiki/Chebyshev_polynomials">Chebyshev 094 * polynomials of the first kind</a> are orthogonal polynomials. 095 * They can be defined by the following recurrence relations:</p><p> 096 * \( 097 * T_0(x) = 1 \\ 098 * T_1(x) = x \\ 099 * T_{k+1}(x) = 2x T_k(x) - T_{k-1}(x) 100 * \) 101 * </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 /** Fixed recurrence coefficients. */ 109 private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE }; 110 /** {@inheritDoc} */ 111 public BigFraction[] generate(int k) { 112 return coeffs; 113 } 114 }); 115 } 116 117 /** 118 * Create a Hermite polynomial. 119 * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite 120 * polynomials</a> are orthogonal polynomials. 121 * They can be defined by the following recurrence relations:</p><p> 122 * \( 123 * H_0(x) = 1 \\ 124 * H_1(x) = 2x \\ 125 * H_{k+1}(x) = 2x H_k(X) - 2k H_{k-1}(x) 126 * \) 127 * </p> 128 129 * @param degree degree of the polynomial 130 * @return Hermite polynomial of specified degree 131 */ 132 public static PolynomialFunction createHermitePolynomial(final int degree) { 133 return buildPolynomial(degree, HERMITE_COEFFICIENTS, 134 new RecurrenceCoefficientsGenerator() { 135 /** {@inheritDoc} */ 136 public BigFraction[] generate(int k) { 137 return new BigFraction[] { 138 BigFraction.ZERO, 139 BigFraction.TWO, 140 new BigFraction(2 * k)}; 141 } 142 }); 143 } 144 145 /** 146 * Create a Laguerre polynomial. 147 * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre 148 * polynomials</a> are orthogonal polynomials. 149 * They can be defined by the following recurrence relations:</p><p> 150 * \( 151 * L_0(x) = 1 \\ 152 * L_1(x) = 1 - x \\ 153 * (k+1) L_{k+1}(x) = (2k + 1 - x) L_k(x) - k L_{k-1}(x) 154 * \) 155 * </p> 156 * @param degree degree of the polynomial 157 * @return Laguerre polynomial of specified degree 158 */ 159 public static PolynomialFunction createLaguerrePolynomial(final int degree) { 160 return buildPolynomial(degree, LAGUERRE_COEFFICIENTS, 161 new RecurrenceCoefficientsGenerator() { 162 /** {@inheritDoc} */ 163 public BigFraction[] generate(int k) { 164 final int kP1 = k + 1; 165 return new BigFraction[] { 166 new BigFraction(2 * k + 1, kP1), 167 new BigFraction(-1, kP1), 168 new BigFraction(k, kP1)}; 169 } 170 }); 171 } 172 173 /** 174 * Create a Legendre polynomial. 175 * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre 176 * polynomials</a> are orthogonal polynomials. 177 * They can be defined by the following recurrence relations:</p><p> 178 * \( 179 * P_0(x) = 1 \\ 180 * P_1(x) = x \\ 181 * (k+1) P_{k+1}(x) = (2k+1) x P_k(x) - k P_{k-1}(x) 182 * \) 183 * </p> 184 * @param degree degree of the polynomial 185 * @return Legendre polynomial of specified degree 186 */ 187 public static PolynomialFunction createLegendrePolynomial(final int degree) { 188 return buildPolynomial(degree, LEGENDRE_COEFFICIENTS, 189 new RecurrenceCoefficientsGenerator() { 190 /** {@inheritDoc} */ 191 public BigFraction[] generate(int k) { 192 final int kP1 = k + 1; 193 return new BigFraction[] { 194 BigFraction.ZERO, 195 new BigFraction(k + kP1, kP1), 196 new BigFraction(k, kP1)}; 197 } 198 }); 199 } 200 201 /** 202 * Create a Jacobi polynomial. 203 * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi 204 * polynomials</a> are orthogonal polynomials. 205 * They can be defined by the following recurrence relations:</p><p> 206 * \( 207 * P_0^{vw}(x) = 1 \\ 208 * P_{-1}^{vw}(x) = 0 \\ 209 * 2k(k + v + w)(2k + v + w - 2) P_k^{vw}(x) = \\ 210 * (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) x + v^2 - w^2] P_{k-1}^{vw}(x) \\ 211 * - 2(k + v - 1)(k + w - 1)(2k + v + w) P_{k-2}^{vw}(x) 212 * \) 213 * </p> 214 * @param degree degree of the polynomial 215 * @param v first exponent 216 * @param w second exponent 217 * @return Jacobi polynomial of specified degree 218 */ 219 public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) { 220 221 // select the appropriate list 222 final JacobiKey key = new JacobiKey(v, w); 223 224 if (!JACOBI_COEFFICIENTS.containsKey(key)) { 225 226 // allocate a new list for v, w 227 final List<BigFraction> list = new ArrayList<BigFraction>(); 228 JACOBI_COEFFICIENTS.put(key, list); 229 230 // Pv,w,0(x) = 1; 231 list.add(BigFraction.ONE); 232 233 // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2 234 list.add(new BigFraction(v - w, 2)); 235 list.add(new BigFraction(2 + v + w, 2)); 236 237 } 238 239 return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key), 240 new RecurrenceCoefficientsGenerator() { 241 /** {@inheritDoc} */ 242 public BigFraction[] generate(int k) { 243 k++; 244 final int kvw = k + v + w; 245 final int twoKvw = kvw + k; 246 final int twoKvwM1 = twoKvw - 1; 247 final int twoKvwM2 = twoKvw - 2; 248 final int den = 2 * k * kvw * twoKvwM2; 249 250 return new BigFraction[] { 251 new BigFraction(twoKvwM1 * (v * v - w * w), den), 252 new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den), 253 new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den) 254 }; 255 } 256 }); 257 258 } 259 260 /** Inner class for Jacobi polynomials keys. */ 261 private static class JacobiKey { 262 263 /** First exponent. */ 264 private final int v; 265 266 /** Second exponent. */ 267 private final int w; 268 269 /** Simple constructor. 270 * @param v first exponent 271 * @param w second exponent 272 */ 273 JacobiKey(final int v, final int w) { 274 this.v = v; 275 this.w = w; 276 } 277 278 /** Get hash code. 279 * @return hash code 280 */ 281 @Override 282 public int hashCode() { 283 return (v << 16) ^ w; 284 } 285 286 /** Check if the instance represent the same key as another instance. 287 * @param key other key 288 * @return true if the instance and the other key refer to the same polynomial 289 */ 290 @Override 291 public boolean equals(final Object key) { 292 293 if ((key == null) || !(key instanceof JacobiKey)) { 294 return false; 295 } 296 297 final JacobiKey otherK = (JacobiKey) key; 298 return (v == otherK.v) && (w == otherK.w); 299 300 } 301 } 302 303 /** 304 * Compute the coefficients of the polynomial \(P_s(x)\) 305 * whose values at point {@code x} will be the same as the those from the 306 * original polynomial \(P(x)\) when computed at {@code x + shift}. 307 * <p> 308 * More precisely, let \(\Delta = \) {@code shift} and let 309 * \(P_s(x) = P(x + \Delta)\). The returned array 310 * consists of the coefficients of \(P_s\). So if \(a_0, ..., a_{n-1}\) 311 * are the coefficients of \(P\), then the returned array 312 * \(b_0, ..., b_{n-1}\) satisfies the identity 313 * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\). 314 * 315 * @param coefficients Coefficients of the original polynomial. 316 * @param shift Shift value. 317 * @return the coefficients \(b_i\) of the shifted 318 * polynomial. 319 */ 320 public static double[] shift(final double[] coefficients, 321 final double shift) { 322 final int dp1 = coefficients.length; 323 final double[] newCoefficients = new double[dp1]; 324 325 // Pascal triangle. 326 final int[][] coeff = new int[dp1][dp1]; 327 for (int i = 0; i < dp1; i++){ 328 for(int j = 0; j <= i; j++){ 329 coeff[i][j] = (int) CombinatoricsUtils.binomialCoefficient(i, j); 330 } 331 } 332 333 // First polynomial coefficient. 334 for (int i = 0; i < dp1; i++){ 335 newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i); 336 } 337 338 // Superior order. 339 final int d = dp1 - 1; 340 for (int i = 0; i < d; i++) { 341 for (int j = i; j < d; j++){ 342 newCoefficients[i + 1] += coeff[j + 1][j - i] * 343 coefficients[j + 1] * FastMath.pow(shift, j - i); 344 } 345 } 346 347 return newCoefficients; 348 } 349 350 351 /** Get the coefficients array for a given degree. 352 * @param degree degree of the polynomial 353 * @param coefficients list where the computed coefficients are stored 354 * @param generator recurrence coefficients generator 355 * @return coefficients array 356 */ 357 private static PolynomialFunction buildPolynomial(final int degree, 358 final List<BigFraction> coefficients, 359 final RecurrenceCoefficientsGenerator generator) { 360 synchronized (coefficients) { 361 final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1; 362 if (degree > maxDegree) { 363 computeUpToDegree(degree, maxDegree, generator, coefficients); 364 } 365 } 366 367 // coefficient for polynomial 0 is l [0] 368 // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1) 369 // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2) 370 // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3) 371 // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4) 372 // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5) 373 // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6) 374 // ... 375 final int start = degree * (degree + 1) / 2; 376 377 final double[] a = new double[degree + 1]; 378 for (int i = 0; i <= degree; ++i) { 379 a[i] = coefficients.get(start + i).doubleValue(); 380 } 381 382 // build the polynomial 383 return new PolynomialFunction(a); 384 385 } 386 387 /** Compute polynomial coefficients up to a given degree. 388 * @param degree maximal degree 389 * @param maxDegree current maximal degree 390 * @param generator recurrence coefficients generator 391 * @param coefficients list where the computed coefficients should be appended 392 */ 393 private static void computeUpToDegree(final int degree, final int maxDegree, 394 final RecurrenceCoefficientsGenerator generator, 395 final List<BigFraction> coefficients) { 396 397 int startK = (maxDegree - 1) * maxDegree / 2; 398 for (int k = maxDegree; k < degree; ++k) { 399 400 // start indices of two previous polynomials Pk(X) and Pk-1(X) 401 int startKm1 = startK; 402 startK += k; 403 404 // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X) 405 BigFraction[] ai = generator.generate(k); 406 407 BigFraction ck = coefficients.get(startK); 408 BigFraction ckm1 = coefficients.get(startKm1); 409 410 // degree 0 coefficient 411 coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2]))); 412 413 // degree 1 to degree k-1 coefficients 414 for (int i = 1; i < k; ++i) { 415 final BigFraction ckPrev = ck; 416 ck = coefficients.get(startK + i); 417 ckm1 = coefficients.get(startKm1 + i); 418 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2]))); 419 } 420 421 // degree k coefficient 422 final BigFraction ckPrev = ck; 423 ck = coefficients.get(startK + k); 424 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1]))); 425 426 // degree k+1 coefficient 427 coefficients.add(ck.multiply(ai[1])); 428 429 } 430 431 } 432 433 /** Interface for recurrence coefficients generation. */ 434 private interface RecurrenceCoefficientsGenerator { 435 /** 436 * Generate recurrence coefficients. 437 * @param k highest degree of the polynomials used in the recurrence 438 * @return an array of three coefficients such that 439 * \( P_{k+1}(x) = (a[0] + a[1] x) P_k(x) - a[2] P_{k-1}(x) \) 440 */ 441 BigFraction[] generate(int k); 442 } 443 444}