PolynomialsUtils.java

  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.  */

  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;

  24. import java.util.ArrayList;
  25. import java.util.HashMap;
  26. import java.util.List;
  27. import java.util.Map;

  29. import org.hipparchus.fraction.BigFraction;
  30. import org.hipparchus.util.CombinatoricsUtils;
  31. import org.hipparchus.util.FastMath;

  33. /**
  34.  * A collection of static methods that operate on or return polynomials.
  35.  *
  36.  */
  37. public class PolynomialsUtils {

  39.     /** Coefficients for Chebyshev polynomials. */
  40.     private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS;

  42.     /** Coefficients for Hermite polynomials. */
  43.     private static final List<BigFraction> HERMITE_COEFFICIENTS;

  45.     /** Coefficients for Laguerre polynomials. */
  46.     private static final List<BigFraction> LAGUERRE_COEFFICIENTS;

  48.     /** Coefficients for Legendre polynomials. */
  49.     private static final List<BigFraction> LEGENDRE_COEFFICIENTS;

  51.     /** Coefficients for Jacobi polynomials. */
  52.     private static final Map<JacobiKey, List<BigFraction>> JACOBI_COEFFICIENTS;

  54.     static {

  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);

  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);

  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);

  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);

  84.         // initialize map for Jacobi polynomials
  85.         JACOBI_COEFFICIENTS = new HashMap<>();

  87.     }

  89.     /**
  90.      * Private constructor, to prevent instantiation.
  91.      */
  92.     private PolynomialsUtils() {
  93.     }

  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.clone();
  118.             }
  119.         });
  120.     }

  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>

  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.     }

  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.     }

  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.     }

  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) {

  229.         // select the appropriate list
  230.         final JacobiKey key = new JacobiKey(v, w);

  232.         if (!JACOBI_COEFFICIENTS.containsKey(key)) {

  234.             // allocate a new list for v, w
  235.             final List<BigFraction> list = new ArrayList<>();
  236.             JACOBI_COEFFICIENTS.put(key, list);

  238.             // Pv,w,0(x) = 1;
  239.             list.add(BigFraction.ONE);

  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));

  245.         }

  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;

  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.         });

  267.     }

  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];

  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.         }

  299.         // First polynomial coefficient.
  300.         for (int i = 0; i < dp1; i++){
  301.             newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i);
  302.         }

  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.         }

  313.         return newCoefficients;
  314.     }


  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) {

  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.         }

  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;

  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.         }

  349.         // build the polynomial
  350.         return new PolynomialFunction(a);

  352.     }

  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) {

  364.         int startK = (maxDegree - 1) * maxDegree / 2;
  365.         for (int k = maxDegree; k < degree; ++k) {

  367.             // start indices of two previous polynomials Pk(X) and Pk-1(X)
  368.             int startKm1 = startK;
  369.             startK += k;

  371.             // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
  372.             BigFraction[] ai = generator.generate(k);

  374.             BigFraction ck     = coefficients.get(startK);
  375.             BigFraction ckm1   = coefficients.get(startKm1);

  377.             // degree 0 coefficient
  378.             coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));

  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.             }

  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])));

  393.             // degree k+1 coefficient
  394.             coefficients.add(ck.multiply(ai[1]));

  396.         }

  398.     }

  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.     }

  411. }
Created with JaCoCo 0.8.12.202403310830