1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math4.legacy.analysis.polynomials;
18
19 import java.util.ArrayList;
20 import java.util.HashMap;
21 import java.util.List;
22 import java.util.Map;
23
24 import org.apache.commons.numbers.fraction.BigFraction;
25 import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
26 import org.apache.commons.math4.core.jdkmath.JdkMath;
27
28 /**
29 * A collection of static methods that operate on or return polynomials.
30 *
31 * @since 2.0
32 */
33 public final class PolynomialsUtils {
34 /** -1. */
35 private static final BigFraction BF_MINUS_ONE = BigFraction.of(-1);
36 /** 2. */
37 private static final BigFraction BF_TWO = BigFraction.of(2);
38
39 /** Coefficients for Chebyshev polynomials. */
40 private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS;
41
42 /** Coefficients for Hermite polynomials. */
43 private static final List<BigFraction> HERMITE_COEFFICIENTS;
44
45 /** Coefficients for Laguerre polynomials. */
46 private static final List<BigFraction> LAGUERRE_COEFFICIENTS;
47
48 /** Coefficients for Legendre polynomials. */
49 private static final List<BigFraction> LEGENDRE_COEFFICIENTS;
50
51 /** Coefficients for Jacobi polynomials. */
52 private static final Map<JacobiKey, List<BigFraction>> JACOBI_COEFFICIENTS;
53
54 static {
55
56 // initialize recurrence for Chebyshev polynomials
57 // T0(X) = 1, T1(X) = 0 + 1 * X
58 CHEBYSHEV_COEFFICIENTS = new ArrayList<>();
59 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
60 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
61 CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
62
63 // initialize recurrence for Hermite polynomials
64 // H0(X) = 1, H1(X) = 0 + 2 * X
65 HERMITE_COEFFICIENTS = new ArrayList<>();
66 HERMITE_COEFFICIENTS.add(BigFraction.ONE);
67 HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
68 HERMITE_COEFFICIENTS.add(BF_TWO);
69
70 // initialize recurrence for Laguerre polynomials
71 // L0(X) = 1, L1(X) = 1 - 1 * X
72 LAGUERRE_COEFFICIENTS = new ArrayList<>();
73 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
74 LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
75 LAGUERRE_COEFFICIENTS.add(BF_MINUS_ONE);
76
77 // initialize recurrence for Legendre polynomials
78 // P0(X) = 1, P1(X) = 0 + 1 * X
79 LEGENDRE_COEFFICIENTS = new ArrayList<>();
80 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
81 LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
82 LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
83
84 // initialize map for Jacobi polynomials
85 JACOBI_COEFFICIENTS = new HashMap<>();
86 }
87
88 /**
89 * Private constructor, to prevent instantiation.
90 */
91 private PolynomialsUtils() {
92 }
93
94 /**
95 * Create a Chebyshev polynomial of the first kind.
96 * <p><a href="https://en.wikipedia.org/wiki/Chebyshev_polynomials">Chebyshev
97 * polynomials of the first kind</a> are orthogonal polynomials.
98 * They can be defined by the following recurrence relations:</p><p>
99 * \(
100 * T_0(x) = 1 \\
101 * T_1(x) = x \\
102 * T_{k+1}(x) = 2x T_k(x) - T_{k-1}(x)
103 * \)
104 * </p>
105 * @param degree degree of the polynomial
106 * @return Chebyshev polynomial of specified degree
107 */
108 public static PolynomialFunction createChebyshevPolynomial(final int degree) {
109 return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
110 new RecurrenceCoefficientsGenerator() {
111 /** Fixed recurrence coefficients. */
112 private final BigFraction[] coeffs = { BigFraction.ZERO, BF_TWO, BigFraction.ONE };
113 /** {@inheritDoc} */
114 @Override
115 public BigFraction[] generate(int k) {
116 return coeffs;
117 }
118 });
119 }
120
121 /**
122 * Create a Hermite polynomial.
123 * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
124 * polynomials</a> are orthogonal polynomials.
125 * They can be defined by the following recurrence relations:</p><p>
126 * \(
127 * H_0(x) = 1 \\
128 * H_1(x) = 2x \\
129 * H_{k+1}(x) = 2x H_k(X) - 2k H_{k-1}(x)
130 * \)
131 * </p>
132
133 * @param degree degree of the polynomial
134 * @return Hermite polynomial of specified degree
135 */
136 public static PolynomialFunction createHermitePolynomial(final int degree) {
137 return buildPolynomial(degree, HERMITE_COEFFICIENTS,
138 new RecurrenceCoefficientsGenerator() {
139 /** {@inheritDoc} */
140 @Override
141 public BigFraction[] generate(int k) {
142 return new BigFraction[] {
143 BigFraction.ZERO,
144 BF_TWO,
145 BigFraction.of(2 * k)};
146 }
147 });
148 }
149
150 /**
151 * Create a Laguerre polynomial.
152 * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
153 * polynomials</a> are orthogonal polynomials.
154 * They can be defined by the following recurrence relations:</p><p>
155 * \(
156 * L_0(x) = 1 \\
157 * L_1(x) = 1 - x \\
158 * (k+1) L_{k+1}(x) = (2k + 1 - x) L_k(x) - k L_{k-1}(x)
159 * \)
160 * </p>
161 * @param degree degree of the polynomial
162 * @return Laguerre polynomial of specified degree
163 */
164 public static PolynomialFunction createLaguerrePolynomial(final int degree) {
165 return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
166 new RecurrenceCoefficientsGenerator() {
167 /** {@inheritDoc} */
168 @Override
169 public BigFraction[] generate(int k) {
170 final int kP1 = k + 1;
171 return new BigFraction[] {
172 BigFraction.of(2 * k + 1, kP1),
173 BigFraction.of(-1, kP1),
174 BigFraction.of(k, kP1)};
175 }
176 });
177 }
178
179 /**
180 * Create a Legendre polynomial.
181 * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
182 * polynomials</a> are orthogonal polynomials.
183 * They can be defined by the following recurrence relations:</p><p>
184 * \(
185 * P_0(x) = 1 \\
186 * P_1(x) = x \\
187 * (k+1) P_{k+1}(x) = (2k+1) x P_k(x) - k P_{k-1}(x)
188 * \)
189 * </p>
190 * @param degree degree of the polynomial
191 * @return Legendre polynomial of specified degree
192 */
193 public static PolynomialFunction createLegendrePolynomial(final int degree) {
194 return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
195 new RecurrenceCoefficientsGenerator() {
196 /** {@inheritDoc} */
197 @Override
198 public BigFraction[] generate(int k) {
199 final int kP1 = k + 1;
200 return new BigFraction[] {
201 BigFraction.ZERO,
202 BigFraction.of(k + kP1, kP1),
203 BigFraction.of(k, kP1)};
204 }
205 });
206 }
207
208 /**
209 * Create a Jacobi polynomial.
210 * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi
211 * polynomials</a> are orthogonal polynomials.
212 * They can be defined by the following recurrence relations:</p><p>
213 * \(
214 * P_0^{vw}(x) = 1 \\
215 * P_{-1}^{vw}(x) = 0 \\
216 * 2k(k + v + w)(2k + v + w - 2) P_k^{vw}(x) = \\
217 * (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) x + v^2 - w^2] P_{k-1}^{vw}(x) \\
218 * - 2(k + v - 1)(k + w - 1)(2k + v + w) P_{k-2}^{vw}(x)
219 * \)
220 * </p>
221 * @param degree degree of the polynomial
222 * @param v first exponent
223 * @param w second exponent
224 * @return Jacobi polynomial of specified degree
225 */
226 public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) {
227
228 // select the appropriate list
229 final JacobiKey key = new JacobiKey(v, w);
230
231 if (!JACOBI_COEFFICIENTS.containsKey(key)) {
232
233 // allocate a new list for v, w
234 final List<BigFraction> list = new ArrayList<>();
235 JACOBI_COEFFICIENTS.put(key, list);
236
237 // Pv,w,0(x) = 1;
238 list.add(BigFraction.ONE);
239
240 // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2
241 list.add(BigFraction.of(v - w, 2));
242 list.add(BigFraction.of(2 + v + w, 2));
243 }
244
245 return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key),
246 new RecurrenceCoefficientsGenerator() {
247 /** {@inheritDoc} */
248 @Override
249 public BigFraction[] generate(int k) {
250 k++;
251 final int kvw = k + v + w;
252 final int twoKvw = kvw + k;
253 final int twoKvwM1 = twoKvw - 1;
254 final int twoKvwM2 = twoKvw - 2;
255 final int den = 2 * k * kvw * twoKvwM2;
256
257 return new BigFraction[] {
258 BigFraction.of(twoKvwM1 * (v * v - w * w), den),
259 BigFraction.of(twoKvwM1 * twoKvw * twoKvwM2, den),
260 BigFraction.of(2 * (k + v - 1) * (k + w - 1) * twoKvw, den)
261 };
262 }
263 });
264 }
265
266 /** Inner class for Jacobi polynomials keys. */
267 private static final class JacobiKey {
268
269 /** First exponent. */
270 private final int v;
271
272 /** Second exponent. */
273 private final int w;
274
275 /** Simple constructor.
276 * @param v first exponent
277 * @param w second exponent
278 */
279 JacobiKey(final int v, final int w) {
280 this.v = v;
281 this.w = w;
282 }
283
284 /** Get hash code.
285 * @return hash code
286 */
287 @Override
288 public int hashCode() {
289 return (v << 16) ^ w;
290 }
291
292 /** Check if the instance represent the same key as another instance.
293 * @param key other key
294 * @return true if the instance and the other key refer to the same polynomial
295 */
296 @Override
297 public boolean equals(final Object key) {
298
299 if (!(key instanceof JacobiKey)) {
300 return false;
301 }
302
303 final JacobiKey otherK = (JacobiKey) key;
304 return v == otherK.v && w == otherK.w;
305 }
306 }
307
308 /**
309 * Compute the coefficients of the polynomial \(P_s(x)\)
310 * whose values at point {@code x} will be the same as the those from the
311 * original polynomial \(P(x)\) when computed at {@code x + shift}.
312 * <p>
313 * More precisely, let \(\Delta = \) {@code shift} and let
314 * \(P_s(x) = P(x + \Delta)\). The returned array
315 * consists of the coefficients of \(P_s\). So if \(a_0, ..., a_{n-1}\)
316 * are the coefficients of \(P\), then the returned array
317 * \(b_0, ..., b_{n-1}\) satisfies the identity
318 * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
319 *
320 * @param coefficients Coefficients of the original polynomial.
321 * @param shift Shift value.
322 * @return the coefficients \(b_i\) of the shifted
323 * polynomial.
324 */
325 public static double[] shift(final double[] coefficients,
326 final double shift) {
327 final int dp1 = coefficients.length;
328 final double[] newCoefficients = new double[dp1];
329
330 // Pascal triangle.
331 final int[][] coeff = new int[dp1][dp1];
332 for (int i = 0; i < dp1; i++){
333 for(int j = 0; j <= i; j++){
334 coeff[i][j] = (int) BinomialCoefficient.value(i, j);
335 }
336 }
337
338 // First polynomial coefficient.
339 for (int i = 0; i < dp1; i++){
340 newCoefficients[0] += coefficients[i] * JdkMath.pow(shift, i);
341 }
342
343 // Superior order.
344 final int d = dp1 - 1;
345 for (int i = 0; i < d; i++) {
346 for (int j = i; j < d; j++){
347 newCoefficients[i + 1] += coeff[j + 1][j - i] *
348 coefficients[j + 1] * JdkMath.pow(shift, j - i);
349 }
350 }
351
352 return newCoefficients;
353 }
354
355
356 /** Get the coefficients array for a given degree.
357 * @param degree degree of the polynomial
358 * @param coefficients list where the computed coefficients are stored
359 * @param generator recurrence coefficients generator
360 * @return coefficients array
361 */
362 private static PolynomialFunction buildPolynomial(final int degree,
363 final List<BigFraction> coefficients,
364 final RecurrenceCoefficientsGenerator generator) {
365 // Synchronizing on a method parameter is not safe; however, in this
366 // case, the lock object is an immutable field that belongs to this
367 // class.
368 synchronized (coefficients) {
369 final int maxDegree = (int) JdkMath.floor(JdkMath.sqrt(2.0 * coefficients.size())) - 1;
370 if (degree > maxDegree) {
371 computeUpToDegree(degree, maxDegree, generator, coefficients);
372 }
373 }
374
375 // coefficient for polynomial 0 is l [0]
376 // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
377 // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
378 // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
379 // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
380 // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
381 // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
382 // ...
383 final int start = degree * (degree + 1) / 2;
384
385 final double[] a = new double[degree + 1];
386 for (int i = 0; i <= degree; ++i) {
387 a[i] = coefficients.get(start + i).doubleValue();
388 }
389
390 // build the polynomial
391 return new PolynomialFunction(a);
392 }
393
394 /** Compute polynomial coefficients up to a given degree.
395 * @param degree maximal degree
396 * @param maxDegree current maximal degree
397 * @param generator recurrence coefficients generator
398 * @param coefficients list where the computed coefficients should be appended
399 */
400 private static void computeUpToDegree(final int degree, final int maxDegree,
401 final RecurrenceCoefficientsGenerator generator,
402 final List<BigFraction> coefficients) {
403
404 int startK = (maxDegree - 1) * maxDegree / 2;
405 for (int k = maxDegree; k < degree; ++k) {
406
407 // start indices of two previous polynomials Pk(X) and Pk-1(X)
408 int startKm1 = startK;
409 startK += k;
410
411 // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
412 BigFraction[] ai = generator.generate(k);
413
414 BigFraction ck = coefficients.get(startK);
415 BigFraction ckm1 = coefficients.get(startKm1);
416
417 // degree 0 coefficient
418 coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
419
420 // degree 1 to degree k-1 coefficients
421 for (int i = 1; i < k; ++i) {
422 final BigFraction ckPrev = ck;
423 ck = coefficients.get(startK + i);
424 ckm1 = coefficients.get(startKm1 + i);
425 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
426 }
427
428 // degree k coefficient
429 final BigFraction ckPrev = ck;
430 ck = coefficients.get(startK + k);
431 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
432
433 // degree k+1 coefficient
434 coefficients.add(ck.multiply(ai[1]));
435 }
436 }
437
438 /** Interface for recurrence coefficients generation. */
439 private interface RecurrenceCoefficientsGenerator {
440 /**
441 * Generate recurrence coefficients.
442 * @param k highest degree of the polynomials used in the recurrence
443 * @return an array of three coefficients such that
444 * \( P_{k+1}(x) = (a[0] + a[1] x) P_k(x) - a[2] P_{k-1}(x) \)
445 */
446 BigFraction[] generate(int k);
447 }
448 }