    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.analysis.polynomials;
    018
    019    import java.util.ArrayList;
    020    import java.util.HashMap;
    021    import java.util.List;
    022    import java.util.Map;
    023
    024    import org.apache.commons.math3.fraction.BigFraction;
    025    import org.apache.commons.math3.util.ArithmeticUtils;
    026    import org.apache.commons.math3.util.FastMath;
    027
    028    /**
    029     * A collection of static methods that operate on or return polynomials.
    030     *
    031     * @version $Id: PolynomialsUtils.java 1364387 2012-07-22 18:14:11Z tn $
    032     * @since 2.0
    033     */
    034    public class PolynomialsUtils {
    035
    036        /** Coefficients for Chebyshev polynomials. */
    037        private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS;
    038
    039        /** Coefficients for Hermite polynomials. */
    040        private static final List<BigFraction> HERMITE_COEFFICIENTS;
    041
    042        /** Coefficients for Laguerre polynomials. */
    043        private static final List<BigFraction> LAGUERRE_COEFFICIENTS;
    044
    045        /** Coefficients for Legendre polynomials. */
    046        private static final List<BigFraction> LEGENDRE_COEFFICIENTS;
    047
    048        /** Coefficients for Jacobi polynomials. */
    049        private static final Map<JacobiKey, List<BigFraction>> JACOBI_COEFFICIENTS;
    050
    051        static {
    052
    053            // initialize recurrence for Chebyshev polynomials
    054            // T0(X) = 1, T1(X) = 0 + 1 * X
    055            CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
    056            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
    057            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
    058            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
    059
    060            // initialize recurrence for Hermite polynomials
    061            // H0(X) = 1, H1(X) = 0 + 2 * X
    062            HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
    063            HERMITE_COEFFICIENTS.add(BigFraction.ONE);
    064            HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
    065            HERMITE_COEFFICIENTS.add(BigFraction.TWO);
    066
    067            // initialize recurrence for Laguerre polynomials
    068            // L0(X) = 1, L1(X) = 1 - 1 * X
    069            LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
    070            LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
    071            LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
    072            LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
    073
    074            // initialize recurrence for Legendre polynomials
    075            // P0(X) = 1, P1(X) = 0 + 1 * X
    076            LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
    077            LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
    078            LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
    079            LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
    080
    081            // initialize map for Jacobi polynomials
    082            JACOBI_COEFFICIENTS = new HashMap<JacobiKey, List<BigFraction>>();
    083
    084        }
    085
    086        /**
    087         * Private constructor, to prevent instantiation.
    088         */
    089        private PolynomialsUtils() {
    090        }
    091
    092        /**
    093         * Create a Chebyshev polynomial of the first kind.
    094         * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
    095         * polynomials of the first kind</a> are orthogonal polynomials.
    096         * They can be defined by the following recurrence relations:
    097         * <pre>
    098         *  T<sub>0</sub>(X)   = 1
    099         *  T<sub>1</sub>(X)   = X
    100         *  T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
    101         * </pre></p>
    102         * @param degree degree of the polynomial
    103         * @return Chebyshev polynomial of specified degree
    104         */
    105        public static PolynomialFunction createChebyshevPolynomial(final int degree) {
    106            return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
    107                    new RecurrenceCoefficientsGenerator() {
    108                private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
    109                /** {@inheritDoc} */
    110                public BigFraction[] generate(int k) {
    111                    return coeffs;
    112                }
    113            });
    114        }
    115
    116        /**
    117         * Create a Hermite polynomial.
    118         * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
    119         * polynomials</a> are orthogonal polynomials.
    120         * They can be defined by the following recurrence relations:
    121         * <pre>
    122         *  H<sub>0</sub>(X)   = 1
    123         *  H<sub>1</sub>(X)   = 2X
    124         *  H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
    125         * </pre></p>
    126
    127         * @param degree degree of the polynomial
    128         * @return Hermite polynomial of specified degree
    129         */
    130        public static PolynomialFunction createHermitePolynomial(final int degree) {
    131            return buildPolynomial(degree, HERMITE_COEFFICIENTS,
    132                    new RecurrenceCoefficientsGenerator() {
    133                /** {@inheritDoc} */
    134                public BigFraction[] generate(int k) {
    135                    return new BigFraction[] {
    136                            BigFraction.ZERO,
    137                            BigFraction.TWO,
    138                            new BigFraction(2 * k)};
    139                }
    140            });
    141        }
    142
    143        /**
    144         * Create a Laguerre polynomial.
    145         * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
    146         * polynomials</a> are orthogonal polynomials.
    147         * They can be defined by the following recurrence relations:
    148         * <pre>
    149         *        L<sub>0</sub>(X)   = 1
    150         *        L<sub>1</sub>(X)   = 1 - X
    151         *  (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
    152         * </pre></p>
    153         * @param degree degree of the polynomial
    154         * @return Laguerre polynomial of specified degree
    155         */
    156        public static PolynomialFunction createLaguerrePolynomial(final int degree) {
    157            return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
    158                    new RecurrenceCoefficientsGenerator() {
    159                /** {@inheritDoc} */
    160                public BigFraction[] generate(int k) {
    161                    final int kP1 = k + 1;
    162                    return new BigFraction[] {
    163                            new BigFraction(2 * k + 1, kP1),
    164                            new BigFraction(-1, kP1),
    165                            new BigFraction(k, kP1)};
    166                }
    167            });
    168        }
    169
    170        /**
    171         * Create a Legendre polynomial.
    172         * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
    173         * polynomials</a> are orthogonal polynomials.
    174         * They can be defined by the following recurrence relations:
    175         * <pre>
    176         *        P<sub>0</sub>(X)   = 1
    177         *        P<sub>1</sub>(X)   = X
    178         *  (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
    179         * </pre></p>
    180         * @param degree degree of the polynomial
    181         * @return Legendre polynomial of specified degree
    182         */
    183        public static PolynomialFunction createLegendrePolynomial(final int degree) {
    184            return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
    185                                   new RecurrenceCoefficientsGenerator() {
    186                /** {@inheritDoc} */
    187                public BigFraction[] generate(int k) {
    188                    final int kP1 = k + 1;
    189                    return new BigFraction[] {
    190                            BigFraction.ZERO,
    191                            new BigFraction(k + kP1, kP1),
    192                            new BigFraction(k, kP1)};
    193                }
    194            });
    195        }
    196
    197        /**
    198         * Create a Jacobi polynomial.
    199         * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi
    200         * polynomials</a> are orthogonal polynomials.
    201         * They can be defined by the following recurrence relations:
    202         * <pre>
    203         *        P<sub>0</sub><sup>vw</sup>(X)   = 1
    204         *        P<sub>-1</sub><sup>vw</sup>(X)  = 0
    205         *  2k(k + v + w)(2k + v + w - 2) P<sub>k</sub><sup>vw</sup>(X) =
    206         *  (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) X + v<sup>2</sup> - w<sup>2</sup>] P<sub>k-1</sub><sup>vw</sup>(X)
    207         *  - 2(k + v - 1)(k + w - 1)(2k + v + w) P<sub>k-2</sub><sup>vw</sup>(X)
    208         * </pre></p>
    209         * @param degree degree of the polynomial
    210         * @param v first exponent
    211         * @param w second exponent
    212         * @return Jacobi polynomial of specified degree
    213         */
    214        public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) {
    215
    216            // select the appropriate list
    217            final JacobiKey key = new JacobiKey(v, w);
    218
    219            if (!JACOBI_COEFFICIENTS.containsKey(key)) {
    220
    221                // allocate a new list for v, w
    222                final List<BigFraction> list = new ArrayList<BigFraction>();
    223                JACOBI_COEFFICIENTS.put(key, list);
    224
    225                // Pv,w,0(x) = 1;
    226                list.add(BigFraction.ONE);
    227
    228                // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2
    229                list.add(new BigFraction(v - w, 2));
    230                list.add(new BigFraction(2 + v + w, 2));
    231
    232            }
    233
    234            return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key),
    235                                   new RecurrenceCoefficientsGenerator() {
    236                /** {@inheritDoc} */
    237                public BigFraction[] generate(int k) {
    238                    k++;
    239                    final int kvw      = k + v + w;
    240                    final int twoKvw   = kvw + k;
    241                    final int twoKvwM1 = twoKvw - 1;
    242                    final int twoKvwM2 = twoKvw - 2;
    243                    final int den      = 2 * k *  kvw * twoKvwM2;
    244
    245                    return new BigFraction[] {
    246                        new BigFraction(twoKvwM1 * (v * v - w * w), den),
    247                        new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den),
    248                        new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den)
    249                    };
    250                }
    251            });
    252
    253        }
    254
    255        /** Inner class for Jacobi polynomials keys. */
    256        private static class JacobiKey {
    257
    258            /** First exponent. */
    259            private final int v;
    260
    261            /** Second exponent. */
    262            private final int w;
    263
    264            /** Simple constructor.
    265             * @param v first exponent
    266             * @param w second exponent
    267             */
    268            public JacobiKey(final int v, final int w) {
    269                this.v = v;
    270                this.w = w;
    271            }
    272
    273            /** Get hash code.
    274             * @return hash code
    275             */
    276            @Override
    277            public int hashCode() {
    278                return (v << 16) ^ w;
    279            }
    280
    281            /** Check if the instance represent the same key as another instance.
    282             * @param key other key
    283             * @return true if the instance and the other key refer to the same polynomial
    284             */
    285            @Override
    286            public boolean equals(final Object key) {
    287
    288                if ((key == null) || !(key instanceof JacobiKey)) {
    289                    return false;
    290                }
    291
    292                final JacobiKey otherK = (JacobiKey) key;
    293                return (v == otherK.v) && (w == otherK.w);
    294
    295            }
    296        }
    297
    298        /**
    299         * Compute the coefficients of the polynomial <code>P<sub>s</sub>(x)</code>
    300         * whose values at point {@code x} will be the same as the those from the
    301         * original polynomial <code>P(x)</code> when computed at {@code x + shift}.
    302         * Thus, if <code>P(x) = &Sigma;<sub>i</sub> a<sub>i</sub> x<sup>i</sup></code>,
    303         * then
    304         * <pre>
    305         *  <table>
    306         *   <tr>
    307         *    <td><code>P<sub>s</sub>(x)</td>
    308         *    <td>= &Sigma;<sub>i</sub> b<sub>i</sub> x<sup>i</sup></code></td>
    309         *   </tr>
    310         *   <tr>
    311         *    <td></td>
    312         *    <td>= &Sigma;<sub>i</sub> a<sub>i</sub> (x + shift)<sup>i</sup></code></td>
    313         *   </tr>
    314         *  </table>
    315         * </pre>
    316         *
    317         * @param coefficients Coefficients of the original polynomial.
    318         * @param shift Shift value.
    319         * @return the coefficients <code>b<sub>i</sub></code> of the shifted
    320         * polynomial.
    321         */
    322        public static double[] shift(final double[] coefficients,
    323                                     final double shift) {
    324            final int dp1 = coefficients.length;
    325            final double[] newCoefficients = new double[dp1];
    326
    327            // Pascal triangle.
    328            final int[][] coeff = new int[dp1][dp1];
    329            for (int i = 0; i < dp1; i++){
    330                for(int j = 0; j <= i; j++){
    331                    coeff[i][j] = (int) ArithmeticUtils.binomialCoefficient(i, j);
    332                }
    333            }
    334
    335            // First polynomial coefficient.
    336            for (int i = 0; i < dp1; i++){
    337                newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i);
    338            }
    339
    340            // Superior order.
    341            final int d = dp1 - 1;
    342            for (int i = 0; i < d; i++) {
    343                for (int j = i; j < d; j++){
    344                    newCoefficients[i + 1] += coeff[j + 1][j - i] *
    345                        coefficients[j + 1] * FastMath.pow(shift, j - i);
    346                }
    347            }
    348
    349            return newCoefficients;
    350        }
    351
    352
    353        /** Get the coefficients array for a given degree.
    354         * @param degree degree of the polynomial
    355         * @param coefficients list where the computed coefficients are stored
    356         * @param generator recurrence coefficients generator
    357         * @return coefficients array
    358         */
    359        private static PolynomialFunction buildPolynomial(final int degree,
    360                                                          final List<BigFraction> coefficients,
    361                                                          final RecurrenceCoefficientsGenerator generator) {
    362
    363            final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
    364            synchronized (PolynomialsUtils.class) {
    365                if (degree > maxDegree) {
    366                    computeUpToDegree(degree, maxDegree, generator, coefficients);
    367                }
    368            }
    369
    370            // coefficient  for polynomial 0 is  l [0]
    371            // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
    372            // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
    373            // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
    374            // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
    375            // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
    376            // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
    377            // ...
    378            final int start = degree * (degree + 1) / 2;
    379
    380            final double[] a = new double[degree + 1];
    381            for (int i = 0; i <= degree; ++i) {
    382                a[i] = coefficients.get(start + i).doubleValue();
    383            }
    384
    385            // build the polynomial
    386            return new PolynomialFunction(a);
    387
    388        }
    389
    390        /** Compute polynomial coefficients up to a given degree.
    391         * @param degree maximal degree
    392         * @param maxDegree current maximal degree
    393         * @param generator recurrence coefficients generator
    394         * @param coefficients list where the computed coefficients should be appended
    395         */
    396        private static void computeUpToDegree(final int degree, final int maxDegree,
    397                                              final RecurrenceCoefficientsGenerator generator,
    398                                              final List<BigFraction> coefficients) {
    399
    400            int startK = (maxDegree - 1) * maxDegree / 2;
    401            for (int k = maxDegree; k < degree; ++k) {
    402
    403                // start indices of two previous polynomials Pk(X) and Pk-1(X)
    404                int startKm1 = startK;
    405                startK += k;
    406
    407                // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
    408                BigFraction[] ai = generator.generate(k);
    409
    410                BigFraction ck     = coefficients.get(startK);
    411                BigFraction ckm1   = coefficients.get(startKm1);
    412
    413                // degree 0 coefficient
    414                coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
    415
    416                // degree 1 to degree k-1 coefficients
    417                for (int i = 1; i < k; ++i) {
    418                    final BigFraction ckPrev = ck;
    419                    ck     = coefficients.get(startK + i);
    420                    ckm1   = coefficients.get(startKm1 + i);
    421                    coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
    422                }
    423
    424                // degree k coefficient
    425                final BigFraction ckPrev = ck;
    426                ck = coefficients.get(startK + k);
    427                coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
    428
    429                // degree k+1 coefficient
    430                coefficients.add(ck.multiply(ai[1]));
    431
    432            }
    433
    434        }
    435
    436        /** Interface for recurrence coefficients generation. */
    437        private interface RecurrenceCoefficientsGenerator {
    438            /**
    439             * Generate recurrence coefficients.
    440             * @param k highest degree of the polynomials used in the recurrence
    441             * @return an array of three coefficients such that
    442             * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
    443             */
    444            BigFraction[] generate(int k);
    445        }
    446
    447    }