Home | History | Annotate | Download | only in polynomials 
      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.math.analysis.polynomials;
     18 
     19 import java.util.ArrayList;
     20 
     21 import org.apache.commons.math.fraction.BigFraction;
     22 import org.apache.commons.math.util.FastMath;
     23 
     24 /**
     25  * A collection of static methods that operate on or return polynomials.
     26  *
     27  * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $
     28  * @since 2.0
     29  */
     30 public class PolynomialsUtils {
     31 
     32     /** Coefficients for Chebyshev polynomials. */
     33     private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
     34 
     35     /** Coefficients for Hermite polynomials. */
     36     private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
     37 
     38     /** Coefficients for Laguerre polynomials. */
     39     private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
     40 
     41     /** Coefficients for Legendre polynomials. */
     42     private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
     43 
     44     static {
     45 
     46         // initialize recurrence for Chebyshev polynomials
     47         // T0(X) = 1, T1(X) = 0 + 1 * X
     48         CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
     49         CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
     50         CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
     51         CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
     52 
     53         // initialize recurrence for Hermite polynomials
     54         // H0(X) = 1, H1(X) = 0 + 2 * X
     55         HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
     56         HERMITE_COEFFICIENTS.add(BigFraction.ONE);
     57         HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
     58         HERMITE_COEFFICIENTS.add(BigFraction.TWO);
     59 
     60         // initialize recurrence for Laguerre polynomials
     61         // L0(X) = 1, L1(X) = 1 - 1 * X
     62         LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
     63         LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
     64         LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
     65         LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
     66 
     67         // initialize recurrence for Legendre polynomials
     68         // P0(X) = 1, P1(X) = 0 + 1 * X
     69         LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
     70         LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
     71         LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
     72         LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
     73 
     74     }
     75 
     76     /**
     77      * Private constructor, to prevent instantiation.
     78      */
     79     private PolynomialsUtils() {
     80     }
     81 
     82     /**
     83      * Create a Chebyshev polynomial of the first kind.
     84      * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
     85      * polynomials of the first kind</a> are orthogonal polynomials.
     86      * They can be defined by the following recurrence relations:
     87      * <pre>
     88      *  T<sub>0</sub>(X)   = 1
     89      *  T<sub>1</sub>(X)   = X
     90      *  T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
     91      * </pre></p>
     92      * @param degree degree of the polynomial
     93      * @return Chebyshev polynomial of specified degree
     94      */
     95     public static PolynomialFunction createChebyshevPolynomial(final int degree) {
     96         return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
     97                 new RecurrenceCoefficientsGenerator() {
     98             private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
     99             /** {@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 }
    282