    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.math.analysis.polynomials;
    018
    019    import java.util.ArrayList;
    020
    021    import org.apache.commons.math.fraction.BigFraction;
    022    import org.apache.commons.math.util.FastMath;
    023
    024    /**
    025     * A collection of static methods that operate on or return polynomials.
    026     *
    027     * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $
    028     * @since 2.0
    029     */
    030    public class PolynomialsUtils {
    031
    032        /** Coefficients for Chebyshev polynomials. */
    033        private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
    034
    035        /** Coefficients for Hermite polynomials. */
    036        private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
    037
    038        /** Coefficients for Laguerre polynomials. */
    039        private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
    040
    041        /** Coefficients for Legendre polynomials. */
    042        private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
    043
    044        static {
    045
    046            // initialize recurrence for Chebyshev polynomials
    047            // T0(X) = 1, T1(X) = 0 + 1 * X
    048            CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
    049            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
    050            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
    051            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
    052
    053            // initialize recurrence for Hermite polynomials
    054            // H0(X) = 1, H1(X) = 0 + 2 * X
    055            HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
    056            HERMITE_COEFFICIENTS.add(BigFraction.ONE);
    057            HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
    058            HERMITE_COEFFICIENTS.add(BigFraction.TWO);
    059
    060            // initialize recurrence for Laguerre polynomials
    061            // L0(X) = 1, L1(X) = 1 - 1 * X
    062            LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
    063            LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
    064            LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
    065            LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
    066
    067            // initialize recurrence for Legendre polynomials
    068            // P0(X) = 1, P1(X) = 0 + 1 * X
    069            LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
    070            LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
    071            LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
    072            LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
    073
    074        }
    075
    076        /**
    077         * Private constructor, to prevent instantiation.
    078         */
    079        private PolynomialsUtils() {
    080        }
    081
    082        /**
    083         * Create a Chebyshev polynomial of the first kind.
    084         * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
    085         * polynomials of the first kind</a> are orthogonal polynomials.
    086         * They can be defined by the following recurrence relations:
    087         * <pre>
    088         *  T<sub>0</sub>(X)   = 1
    089         *  T<sub>1</sub>(X)   = X
    090         *  T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
    091         * </pre></p>
    092         * @param degree degree of the polynomial
    093         * @return Chebyshev polynomial of specified degree
    094         */
    095        public static PolynomialFunction createChebyshevPolynomial(final int degree) {
    096            return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
    097                    new RecurrenceCoefficientsGenerator() {
    098                private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
    099                /** {@inheritDoc} */
    100                public BigFraction[] generate(int k) {
    101                    return coeffs;
    102                }
    103            });
    104        }
    105
    106        /**
    107         * Create a Hermite polynomial.
    108         * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
    109         * polynomials</a> are orthogonal polynomials.
    110         * They can be defined by the following recurrence relations:
    111         * <pre>
    112         *  H<sub>0</sub>(X)   = 1
    113         *  H<sub>1</sub>(X)   = 2X
    114         *  H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
    115         * </pre></p>
    116
    117         * @param degree degree of the polynomial
    118         * @return Hermite polynomial of specified degree
    119         */
    120        public static PolynomialFunction createHermitePolynomial(final int degree) {
    121            return buildPolynomial(degree, HERMITE_COEFFICIENTS,
    122                    new RecurrenceCoefficientsGenerator() {
    123                /** {@inheritDoc} */
    124                public BigFraction[] generate(int k) {
    125                    return new BigFraction[] {
    126                            BigFraction.ZERO,
    127                            BigFraction.TWO,
    128                            new BigFraction(2 * k)};
    129                }
    130            });
    131        }
    132
    133        /**
    134         * Create a Laguerre polynomial.
    135         * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
    136         * polynomials</a> are orthogonal polynomials.
    137         * They can be defined by the following recurrence relations:
    138         * <pre>
    139         *        L<sub>0</sub>(X)   = 1
    140         *        L<sub>1</sub>(X)   = 1 - X
    141         *  (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
    142         * </pre></p>
    143         * @param degree degree of the polynomial
    144         * @return Laguerre polynomial of specified degree
    145         */
    146        public static PolynomialFunction createLaguerrePolynomial(final int degree) {
    147            return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
    148                    new RecurrenceCoefficientsGenerator() {
    149                /** {@inheritDoc} */
    150                public BigFraction[] generate(int k) {
    151                    final int kP1 = k + 1;
    152                    return new BigFraction[] {
    153                            new BigFraction(2 * k + 1, kP1),
    154                            new BigFraction(-1, kP1),
    155                            new BigFraction(k, kP1)};
    156                }
    157            });
    158        }
    159
    160        /**
    161         * Create a Legendre polynomial.
    162         * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
    163         * polynomials</a> are orthogonal polynomials.
    164         * They can be defined by the following recurrence relations:
    165         * <pre>
    166         *        P<sub>0</sub>(X)   = 1
    167         *        P<sub>1</sub>(X)   = X
    168         *  (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
    169         * </pre></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                public BigFraction[] generate(int k) {
    178                    final int kP1 = k + 1;
    179                    return new BigFraction[] {
    180                            BigFraction.ZERO,
    181                            new BigFraction(k + kP1, kP1),
    182                            new BigFraction(k, kP1)};
    183                }
    184            });
    185        }
    186
    187        /** Get the coefficients array for a given degree.
    188         * @param degree degree of the polynomial
    189         * @param coefficients list where the computed coefficients are stored
    190         * @param generator recurrence coefficients generator
    191         * @return coefficients array
    192         */
    193        private static PolynomialFunction buildPolynomial(final int degree,
    194                                                          final ArrayList<BigFraction> coefficients,
    195                                                          final RecurrenceCoefficientsGenerator generator) {
    196
    197            final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
    198            synchronized (PolynomialsUtils.class) {
    199                if (degree > maxDegree) {
    200                    computeUpToDegree(degree, maxDegree, generator, coefficients);
    201                }
    202            }
    203
    204            // coefficient  for polynomial 0 is  l [0]
    205            // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
    206            // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
    207            // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
    208            // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
    209            // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
    210            // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
    211            // ...
    212            final int start = degree * (degree + 1) / 2;
    213
    214            final double[] a = new double[degree + 1];
    215            for (int i = 0; i <= degree; ++i) {
    216                a[i] = coefficients.get(start + i).doubleValue();
    217            }
    218
    219            // build the polynomial
    220            return new PolynomialFunction(a);
    221
    222        }
    223
    224        /** Compute polynomial coefficients up to a given degree.
    225         * @param degree maximal degree
    226         * @param maxDegree current maximal degree
    227         * @param generator recurrence coefficients generator
    228         * @param coefficients list where the computed coefficients should be appended
    229         */
    230        private static void computeUpToDegree(final int degree, final int maxDegree,
    231                                              final RecurrenceCoefficientsGenerator generator,
    232                                              final ArrayList<BigFraction> coefficients) {
    233
    234            int startK = (maxDegree - 1) * maxDegree / 2;
    235            for (int k = maxDegree; k < degree; ++k) {
    236
    237                // start indices of two previous polynomials Pk(X) and Pk-1(X)
    238                int startKm1 = startK;
    239                startK += k;
    240
    241                // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
    242                BigFraction[] ai = generator.generate(k);
    243
    244                BigFraction ck     = coefficients.get(startK);
    245                BigFraction ckm1   = coefficients.get(startKm1);
    246
    247                // degree 0 coefficient
    248                coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
    249
    250                // degree 1 to degree k-1 coefficients
    251                for (int i = 1; i < k; ++i) {
    252                    final BigFraction ckPrev = ck;
    253                    ck     = coefficients.get(startK + i);
    254                    ckm1   = coefficients.get(startKm1 + i);
    255                    coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
    256                }
    257
    258                // degree k coefficient
    259                final BigFraction ckPrev = ck;
    260                ck = coefficients.get(startK + k);
    261                coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
    262
    263                // degree k+1 coefficient
    264                coefficients.add(ck.multiply(ai[1]));
    265
    266            }
    267
    268        }
    269
    270        /** Interface for recurrence coefficients generation. */
    271        private static interface RecurrenceCoefficientsGenerator {
    272            /**
    273             * Generate recurrence coefficients.
    274             * @param k highest degree of the polynomials used in the recurrence
    275             * @return an array of three coefficients such that
    276             * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
    277             */
    278            BigFraction[] generate(int k);
    279        }
    280
    281    }