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