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.  *      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. import java.util.ArrayList;
  19. import java.util.HashMap;
  20. import java.util.List;
  21. import java.util.Map;

  22. import org.apache.commons.numbers.fraction.BigFraction;
  23. import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
  24. import org.apache.commons.math4.core.jdkmath.JdkMath;

  25. /**
  26.  * A collection of static methods that operate on or return polynomials.
  27.  *
  28.  * @since 2.0
  29.  */
  30. public final class PolynomialsUtils {
  31.     /** -1. */
  32.     private static final BigFraction BF_MINUS_ONE = BigFraction.of(-1);
  33.     /** 2. */
  34.     private static final BigFraction BF_TWO = BigFraction.of(2);

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

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

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

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

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

  45.     static {

  46.         // initialize recurrence for Chebyshev polynomials
  47.         // T0(X) = 1, T1(X) = 0 + 1 * X
  48.         CHEBYSHEV_COEFFICIENTS = new ArrayList<>();
  49.         CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
  50.         CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
  51.         CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);

  52.         // initialize recurrence for Hermite polynomials
  53.         // H0(X) = 1, H1(X) = 0 + 2 * X
  54.         HERMITE_COEFFICIENTS = new ArrayList<>();
  55.         HERMITE_COEFFICIENTS.add(BigFraction.ONE);
  56.         HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
  57.         HERMITE_COEFFICIENTS.add(BF_TWO);

  58.         // initialize recurrence for Laguerre polynomials
  59.         // L0(X) = 1, L1(X) = 1 - 1 * X
  60.         LAGUERRE_COEFFICIENTS = new ArrayList<>();
  61.         LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
  62.         LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
  63.         LAGUERRE_COEFFICIENTS.add(BF_MINUS_ONE);

  64.         // initialize recurrence for Legendre polynomials
  65.         // P0(X) = 1, P1(X) = 0 + 1 * X
  66.         LEGENDRE_COEFFICIENTS = new ArrayList<>();
  67.         LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
  68.         LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
  69.         LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);

  70.         // initialize map for Jacobi polynomials
  71.         JACOBI_COEFFICIENTS = new HashMap<>();
  72.     }

  73.     /**
  74.      * Private constructor, to prevent instantiation.
  75.      */
  76.     private PolynomialsUtils() {
  77.     }

  78.     /**
  79.      * Create a Chebyshev polynomial of the first kind.
  80.      * <p><a href="https://en.wikipedia.org/wiki/Chebyshev_polynomials">Chebyshev
  81.      * polynomials of the first kind</a> are orthogonal polynomials.
  82.      * They can be defined by the following recurrence relations:</p><p>
  83.      * \(
  84.      *    T_0(x) = 1 \\
  85.      *    T_1(x) = x \\
  86.      *    T_{k+1}(x) = 2x T_k(x) - T_{k-1}(x)
  87.      * \)
  88.      * </p>
  89.      * @param degree degree of the polynomial
  90.      * @return Chebyshev polynomial of specified degree
  91.      */
  92.     public static PolynomialFunction createChebyshevPolynomial(final int degree) {
  93.         return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
  94.                 new RecurrenceCoefficientsGenerator() {
  95.             /** Fixed recurrence coefficients. */
  96.             private final BigFraction[] coeffs = { BigFraction.ZERO, BF_TWO, BigFraction.ONE };
  97.             /** {@inheritDoc} */
  98.             @Override
  99.             public BigFraction[] generate(int k) {
  100.                 return coeffs;
  101.             }
  102.         });
  103.     }

  104.     /**
  105.      * Create a Hermite polynomial.
  106.      * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
  107.      * polynomials</a> are orthogonal polynomials.
  108.      * They can be defined by the following recurrence relations:</p><p>
  109.      * \(
  110.      *  H_0(x) = 1 \\
  111.      *  H_1(x) = 2x \\
  112.      *  H_{k+1}(x) = 2x H_k(X) - 2k H_{k-1}(x)
  113.      * \)
  114.      * </p>

  115.      * @param degree degree of the polynomial
  116.      * @return Hermite polynomial of specified degree
  117.      */
  118.     public static PolynomialFunction createHermitePolynomial(final int degree) {
  119.         return buildPolynomial(degree, HERMITE_COEFFICIENTS,
  120.                 new RecurrenceCoefficientsGenerator() {
  121.             /** {@inheritDoc} */
  122.             @Override
  123.             public BigFraction[] generate(int k) {
  124.                 return new BigFraction[] {
  125.                         BigFraction.ZERO,
  126.                         BF_TWO,
  127.                         BigFraction.of(2 * k)};
  128.             }
  129.         });
  130.     }

  131.     /**
  132.      * Create a Laguerre polynomial.
  133.      * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
  134.      * polynomials</a> are orthogonal polynomials.
  135.      * They can be defined by the following recurrence relations:</p><p>
  136.      * \(
  137.      *   L_0(x) = 1 \\
  138.      *   L_1(x) = 1 - x \\
  139.      *   (k+1) L_{k+1}(x) = (2k + 1 - x) L_k(x) - k L_{k-1}(x)
  140.      * \)
  141.      * </p>
  142.      * @param degree degree of the polynomial
  143.      * @return Laguerre polynomial of specified degree
  144.      */
  145.     public static PolynomialFunction createLaguerrePolynomial(final int degree) {
  146.         return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
  147.                 new RecurrenceCoefficientsGenerator() {
  148.             /** {@inheritDoc} */
  149.             @Override
  150.             public BigFraction[] generate(int k) {
  151.                 final int kP1 = k + 1;
  152.                 return new BigFraction[] {
  153.                         BigFraction.of(2 * k + 1, kP1),
  154.                         BigFraction.of(-1, kP1),
  155.                         BigFraction.of(k, kP1)};
  156.             }
  157.         });
  158.     }

  159.     /**
  160.      * Create a Legendre polynomial.
  161.      * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
  162.      * polynomials</a> are orthogonal polynomials.
  163.      * They can be defined by the following recurrence relations:</p><p>
  164.      * \(
  165.      *   P_0(x) = 1 \\
  166.      *   P_1(x) = x \\
  167.      *   (k+1) P_{k+1}(x) = (2k+1) x P_k(x) - k P_{k-1}(x)
  168.      * \)
  169.      * </p>
  170.      * @param degree degree of the polynomial
  171.      * @return Legendre polynomial of specified degree
  172.      */
  173.     public static PolynomialFunction createLegendrePolynomial(final int degree) {
  174.         return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
  175.                                new RecurrenceCoefficientsGenerator() {
  176.             /** {@inheritDoc} */
  177.             @Override
  178.             public BigFraction[] generate(int k) {
  179.                 final int kP1 = k + 1;
  180.                 return new BigFraction[] {
  181.                         BigFraction.ZERO,
  182.                         BigFraction.of(k + kP1, kP1),
  183.                         BigFraction.of(k, kP1)};
  184.             }
  185.         });
  186.     }

  187.     /**
  188.      * Create a Jacobi polynomial.
  189.      * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi
  190.      * polynomials</a> are orthogonal polynomials.
  191.      * They can be defined by the following recurrence relations:</p><p>
  192.      * \(
  193.      *    P_0^{vw}(x) = 1 \\
  194.      *    P_{-1}^{vw}(x) = 0 \\
  195.      *    2k(k + v + w)(2k + v + w - 2) P_k^{vw}(x) = \\
  196.      *    (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) x + v^2 - w^2] P_{k-1}^{vw}(x) \\
  197.      *  - 2(k + v - 1)(k + w - 1)(2k + v + w) P_{k-2}^{vw}(x)
  198.      * \)
  199.      * </p>
  200.      * @param degree degree of the polynomial
  201.      * @param v first exponent
  202.      * @param w second exponent
  203.      * @return Jacobi polynomial of specified degree
  204.      */
  205.     public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) {

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

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

  209.             // allocate a new list for v, w
  210.             final List<BigFraction> list = new ArrayList<>();
  211.             JACOBI_COEFFICIENTS.put(key, list);

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

  214.             // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2
  215.             list.add(BigFraction.of(v - w, 2));
  216.             list.add(BigFraction.of(2 + v + w, 2));
  217.         }

  218.         return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key),
  219.                                new RecurrenceCoefficientsGenerator() {
  220.             /** {@inheritDoc} */
  221.             @Override
  222.             public BigFraction[] generate(int k) {
  223.                 k++;
  224.                 final int kvw      = k + v + w;
  225.                 final int twoKvw   = kvw + k;
  226.                 final int twoKvwM1 = twoKvw - 1;
  227.                 final int twoKvwM2 = twoKvw - 2;
  228.                 final int den      = 2 * k *  kvw * twoKvwM2;

  229.                 return new BigFraction[] {
  230.                     BigFraction.of(twoKvwM1 * (v * v - w * w), den),
  231.                     BigFraction.of(twoKvwM1 * twoKvw * twoKvwM2, den),
  232.                     BigFraction.of(2 * (k + v - 1) * (k + w - 1) * twoKvw, den)
  233.                 };
  234.             }
  235.         });
  236.     }

  237.     /** Inner class for Jacobi polynomials keys. */
  238.     private static final class JacobiKey {

  239.         /** First exponent. */
  240.         private final int v;

  241.         /** Second exponent. */
  242.         private final int w;

  243.         /** Simple constructor.
  244.          * @param v first exponent
  245.          * @param w second exponent
  246.          */
  247.         JacobiKey(final int v, final int w) {
  248.             this.v = v;
  249.             this.w = w;
  250.         }

  251.         /** Get hash code.
  252.          * @return hash code
  253.          */
  254.         @Override
  255.         public int hashCode() {
  256.             return (v << 16) ^ w;
  257.         }

  258.         /** Check if the instance represent the same key as another instance.
  259.          * @param key other key
  260.          * @return true if the instance and the other key refer to the same polynomial
  261.          */
  262.         @Override
  263.         public boolean equals(final Object key) {

  264.             if (!(key instanceof JacobiKey)) {
  265.                 return false;
  266.             }

  267.             final JacobiKey otherK = (JacobiKey) key;
  268.             return v == otherK.v && w == otherK.w;
  269.         }
  270.     }

  271.     /**
  272.      * Compute the coefficients of the polynomial \(P_s(x)\)
  273.      * whose values at point {@code x} will be the same as the those from the
  274.      * original polynomial \(P(x)\) when computed at {@code x + shift}.
  275.      * <p>
  276.      * More precisely, let \(\Delta = \) {@code shift} and let
  277.      * \(P_s(x) = P(x + \Delta)\).  The returned array
  278.      * consists of the coefficients of \(P_s\).  So if \(a_0, ..., a_{n-1}\)
  279.      * are the coefficients of \(P\), then the returned array
  280.      * \(b_0, ..., b_{n-1}\) satisfies the identity
  281.      * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
  282.      *
  283.      * @param coefficients Coefficients of the original polynomial.
  284.      * @param shift Shift value.
  285.      * @return the coefficients \(b_i\) of the shifted
  286.      * polynomial.
  287.      */
  288.     public static double[] shift(final double[] coefficients,
  289.                                  final double shift) {
  290.         final int dp1 = coefficients.length;
  291.         final double[] newCoefficients = new double[dp1];

  292.         // Pascal triangle.
  293.         final int[][] coeff = new int[dp1][dp1];
  294.         for (int i = 0; i < dp1; i++){
  295.             for(int j = 0; j <= i; j++){
  296.                 coeff[i][j] = (int) BinomialCoefficient.value(i, j);
  297.             }
  298.         }

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

  303.         // Superior order.
  304.         final int d = dp1 - 1;
  305.         for (int i = 0; i < d; i++) {
  306.             for (int j = i; j < d; j++){
  307.                 newCoefficients[i + 1] += coeff[j + 1][j - i] *
  308.                     coefficients[j + 1] * JdkMath.pow(shift, j - i);
  309.             }
  310.         }

  311.         return newCoefficients;
  312.     }


  313.     /** Get the coefficients array for a given degree.
  314.      * @param degree degree of the polynomial
  315.      * @param coefficients list where the computed coefficients are stored
  316.      * @param generator recurrence coefficients generator
  317.      * @return coefficients array
  318.      */
  319.     private static PolynomialFunction buildPolynomial(final int degree,
  320.                                                       final List<BigFraction> coefficients,
  321.                                                       final RecurrenceCoefficientsGenerator generator) {
  322.         // Synchronizing on a method parameter is not safe; however, in this
  323.         // case, the lock object is an immutable field that belongs to this
  324.         // class.
  325.         synchronized (coefficients) {
  326.             final int maxDegree = (int) JdkMath.floor(JdkMath.sqrt(2.0 * coefficients.size())) - 1;
  327.             if (degree > maxDegree) {
  328.                 computeUpToDegree(degree, maxDegree, generator, coefficients);
  329.             }
  330.         }

  331.         // coefficient  for polynomial 0 is  l [0]
  332.         // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
  333.         // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
  334.         // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
  335.         // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
  336.         // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
  337.         // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
  338.         // ...
  339.         final int start = degree * (degree + 1) / 2;

  340.         final double[] a = new double[degree + 1];
  341.         for (int i = 0; i <= degree; ++i) {
  342.             a[i] = coefficients.get(start + i).doubleValue();
  343.         }

  344.         // build the polynomial
  345.         return new PolynomialFunction(a);
  346.     }

  347.     /** Compute polynomial coefficients up to a given degree.
  348.      * @param degree maximal degree
  349.      * @param maxDegree current maximal degree
  350.      * @param generator recurrence coefficients generator
  351.      * @param coefficients list where the computed coefficients should be appended
  352.      */
  353.     private static void computeUpToDegree(final int degree, final int maxDegree,
  354.                                           final RecurrenceCoefficientsGenerator generator,
  355.                                           final List<BigFraction> coefficients) {

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

  358.             // start indices of two previous polynomials Pk(X) and Pk-1(X)
  359.             int startKm1 = startK;
  360.             startK += k;

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

  363.             BigFraction ck     = coefficients.get(startK);
  364.             BigFraction ckm1   = coefficients.get(startKm1);

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

  367.             // degree 1 to degree k-1 coefficients
  368.             for (int i = 1; i < k; ++i) {
  369.                 final BigFraction ckPrev = ck;
  370.                 ck     = coefficients.get(startK + i);
  371.                 ckm1   = coefficients.get(startKm1 + i);
  372.                 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
  373.             }

  374.             // degree k coefficient
  375.             final BigFraction ckPrev = ck;
  376.             ck = coefficients.get(startK + k);
  377.             coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));

  378.             // degree k+1 coefficient
  379.             coefficients.add(ck.multiply(ai[1]));
  380.         }
  381.     }

  382.     /** Interface for recurrence coefficients generation. */
  383.     private interface RecurrenceCoefficientsGenerator {
  384.         /**
  385.          * Generate recurrence coefficients.
  386.          * @param k highest degree of the polynomials used in the recurrence
  387.          * @return an array of three coefficients such that
  388.          * \( P_{k+1}(x) = (a[0] + a[1] x) P_k(x) - a[2] P_{k-1}(x) \)
  389.          */
  390.         BigFraction[] generate(int k);
  391.     }
  392. }