View Javadoc
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     /** Inner class for Jacobi polynomials keys. */
270     private static class JacobiKey {
271 
272         /** First exponent. */
273         private final int v;
274 
275         /** Second exponent. */
276         private final int w;
277 
278         /** Simple constructor.
279          * @param v first exponent
280          * @param w second exponent
281          */
282         JacobiKey(final int v, final int w) {
283             this.v = v;
284             this.w = w;
285         }
286 
287         /** Get hash code.
288          * @return hash code
289          */
290         @Override
291         public int hashCode() {
292             return (v << 16) ^ w;
293         }
294 
295         /** Check if the instance represent the same key as another instance.
296          * @param key other key
297          * @return true if the instance and the other key refer to the same polynomial
298          */
299         @Override
300         public boolean equals(final Object key) {
301 
302             if ((key == null) || !(key instanceof JacobiKey)) {
303                 return false;
304             }
305 
306             final JacobiKey otherK = (JacobiKey) key;
307             return (v == otherK.v) && (w == otherK.w);
308 
309         }
310     }
311 
312     /**
313      * Compute the coefficients of the polynomial \(P_s(x)\)
314      * whose values at point {@code x} will be the same as the those from the
315      * original polynomial \(P(x)\) when computed at {@code x + shift}.
316      * <p>
317      * More precisely, let \(\Delta = \) {@code shift} and let
318      * \(P_s(x) = P(x + \Delta)\).  The returned array
319      * consists of the coefficients of \(P_s\).  So if \(a_0, ..., a_{n-1}\)
320      * are the coefficients of \(P\), then the returned array
321      * \(b_0, ..., b_{n-1}\) satisfies the identity
322      * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
323      *
324      * @param coefficients Coefficients of the original polynomial.
325      * @param shift Shift value.
326      * @return the coefficients \(b_i\) of the shifted
327      * polynomial.
328      */
329     public static double[] shift(final double[] coefficients,
330                                  final double shift) {
331         final int dp1 = coefficients.length;
332         final double[] newCoefficients = new double[dp1];
333 
334         // Pascal triangle.
335         final int[][] coeff = new int[dp1][dp1];
336         for (int i = 0; i < dp1; i++){
337             for(int j = 0; j <= i; j++){
338                 coeff[i][j] = (int) CombinatoricsUtils.binomialCoefficient(i, j);
339             }
340         }
341 
342         // First polynomial coefficient.
343         for (int i = 0; i < dp1; i++){
344             newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i);
345         }
346 
347         // Superior order.
348         final int d = dp1 - 1;
349         for (int i = 0; i < d; i++) {
350             for (int j = i; j < d; j++){
351                 newCoefficients[i + 1] += coeff[j + 1][j - i] *
352                     coefficients[j + 1] * FastMath.pow(shift, j - i);
353             }
354         }
355 
356         return newCoefficients;
357     }
358 
359 
360     /** Get the coefficients array for a given degree.
361      * @param degree degree of the polynomial
362      * @param coefficients list where the computed coefficients are stored
363      * @param generator recurrence coefficients generator
364      * @return coefficients array
365      */
366     private static PolynomialFunction buildPolynomial(final int degree,
367                                                       final List<BigFraction> coefficients,
368                                                       final RecurrenceCoefficientsGenerator generator) {
369 
370         synchronized (coefficients) {
371             final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
372             if (degree > maxDegree) {
373                 computeUpToDegree(degree, maxDegree, generator, coefficients);
374             }
375         }
376 
377         // coefficient  for polynomial 0 is  l [0]
378         // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
379         // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
380         // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
381         // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
382         // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
383         // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
384         // ...
385         final int start = degree * (degree + 1) / 2;
386 
387         final double[] a = new double[degree + 1];
388         for (int i = 0; i <= degree; ++i) {
389             a[i] = coefficients.get(start + i).doubleValue();
390         }
391 
392         // build the polynomial
393         return new PolynomialFunction(a);
394 
395     }
396 
397     /** Compute polynomial coefficients up to a given degree.
398      * @param degree maximal degree
399      * @param maxDegree current maximal degree
400      * @param generator recurrence coefficients generator
401      * @param coefficients list where the computed coefficients should be appended
402      */
403     private static void computeUpToDegree(final int degree, final int maxDegree,
404                                           final RecurrenceCoefficientsGenerator generator,
405                                           final List<BigFraction> coefficients) {
406 
407         int startK = (maxDegree - 1) * maxDegree / 2;
408         for (int k = maxDegree; k < degree; ++k) {
409 
410             // start indices of two previous polynomials Pk(X) and Pk-1(X)
411             int startKm1 = startK;
412             startK += k;
413 
414             // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
415             BigFraction[] ai = generator.generate(k);
416 
417             BigFraction ck     = coefficients.get(startK);
418             BigFraction ckm1   = coefficients.get(startKm1);
419 
420             // degree 0 coefficient
421             coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
422 
423             // degree 1 to degree k-1 coefficients
424             for (int i = 1; i < k; ++i) {
425                 final BigFraction ckPrev = ck;
426                 ck     = coefficients.get(startK + i);
427                 ckm1   = coefficients.get(startKm1 + i);
428                 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
429             }
430 
431             // degree k coefficient
432             final BigFraction ckPrev = ck;
433             ck = coefficients.get(startK + k);
434             coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
435 
436             // degree k+1 coefficient
437             coefficients.add(ck.multiply(ai[1]));
438 
439         }
440 
441     }
442 
443     /** Interface for recurrence coefficients generation. */
444     private interface RecurrenceCoefficientsGenerator {
445         /**
446          * Generate recurrence coefficients.
447          * @param k highest degree of the polynomials used in the recurrence
448          * @return an array of three coefficients such that
449          * \( P_{k+1}(x) = (a[0] + a[1] x) P_k(x) - a[2] P_{k-1}(x) \)
450          */
451         BigFraction[] generate(int k);
452     }
453 
454 }