View Javadoc

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  
23  /**
24   * A collection of static methods that operate on or return polynomials.
25   * 
26   * @version $Revision: 760901 $ $Date: 2009-04-01 10:29:18 -0400 (Wed, 01 Apr 2009) $
27   * @since 2.0
28   */
29  public class PolynomialsUtils {
30  
31      /** Coefficients for Chebyshev polynomials. */
32      private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
33  
34      /** Coefficients for Hermite polynomials. */
35      private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
36  
37      /** Coefficients for Laguerre polynomials. */
38      private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
39  
40      /** Coefficients for Legendre polynomials. */
41      private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
42  
43      static {
44  
45          // initialize recurrence for Chebyshev polynomials
46          // T0(X) = 1, T1(X) = 0 + 1 * X
47          CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
48          CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
49          CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
50          CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
51  
52          // initialize recurrence for Hermite polynomials
53          // H0(X) = 1, H1(X) = 0 + 2 * X
54          HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
55          HERMITE_COEFFICIENTS.add(BigFraction.ONE);
56          HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
57          HERMITE_COEFFICIENTS.add(BigFraction.TWO);
58  
59          // initialize recurrence for Laguerre polynomials
60          // L0(X) = 1, L1(X) = 1 - 1 * X
61          LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
62          LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
63          LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
64          LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
65  
66          // initialize recurrence for Legendre polynomials
67          // P0(X) = 1, P1(X) = 0 + 1 * X
68          LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
69          LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
70          LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
71          LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
72  
73      }
74  
75      /**
76       * Private constructor, to prevent instantiation.
77       */
78      private PolynomialsUtils() {
79      }
80  
81      /**
82       * Create a Chebyshev polynomial of the first kind.
83       * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
84       * polynomials of the first kind</a> are orthogonal polynomials.
85       * They can be defined by the following recurrence relations:
86       * <pre>
87       *  T<sub>0</sub>(X)   = 1
88       *  T<sub>1</sub>(X)   = X
89       *  T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
90       * </pre></p>
91       * @param degree degree of the polynomial
92       * @return Chebyshev polynomial of specified degree
93       */
94      public static PolynomialFunction createChebyshevPolynomial(final int degree) {
95          return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
96                  new RecurrenceCoefficientsGenerator() {
97              private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
98              /** {@inheritDoc} */
99              public BigFraction[] generate(int k) {
100                 return coeffs;
101             }
102         });
103     }
104 
105     /**
106      * Create a Hermite polynomial.
107      * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
108      * polynomials</a> are orthogonal polynomials.
109      * They can be defined by the following recurrence relations:
110      * <pre>
111      *  H<sub>0</sub>(X)   = 1
112      *  H<sub>1</sub>(X)   = 2X
113      *  H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
114      * </pre></p>
115 
116      * @param degree degree of the polynomial
117      * @return Hermite polynomial of specified degree
118      */
119     public static PolynomialFunction createHermitePolynomial(final int degree) {
120         return buildPolynomial(degree, HERMITE_COEFFICIENTS,
121                 new RecurrenceCoefficientsGenerator() {
122             /** {@inheritDoc} */
123             public BigFraction[] generate(int k) {
124                 return new BigFraction[] {
125                         BigFraction.ZERO,
126                         BigFraction.TWO,
127                         new BigFraction(2 * k)};
128             }
129         });
130     }
131 
132     /**
133      * Create a Laguerre polynomial.
134      * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
135      * polynomials</a> are orthogonal polynomials.
136      * They can be defined by the following recurrence relations:
137      * <pre>
138      *        L<sub>0</sub>(X)   = 1
139      *        L<sub>1</sub>(X)   = 1 - X
140      *  (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
141      * </pre></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             public BigFraction[] generate(int k) {
150                 final int kP1 = k + 1;
151                 return new BigFraction[] {
152                         new BigFraction(2 * k + 1, kP1),
153                         new BigFraction(-1, kP1),
154                         new BigFraction(k, kP1)};
155             }
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:
164      * <pre>
165      *        P<sub>0</sub>(X)   = 1
166      *        P<sub>1</sub>(X)   = X
167      *  (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
168      * </pre></p>
169      * @param degree degree of the polynomial
170      * @return Legendre polynomial of specified degree
171      */
172     public static PolynomialFunction createLegendrePolynomial(final int degree) {
173         return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
174                                new RecurrenceCoefficientsGenerator() {
175             /** {@inheritDoc} */
176             public BigFraction[] generate(int k) {
177                 final int kP1 = k + 1;
178                 return new BigFraction[] {
179                         BigFraction.ZERO,
180                         new BigFraction(k + kP1, kP1),
181                         new BigFraction(k, kP1)};
182             }
183         });
184     }
185 
186     /** Get the coefficients array for a given degree.
187      * @param degree degree of the polynomial
188      * @param coefficients list where the computed coefficients are stored
189      * @param generator recurrence coefficients generator
190      * @return coefficients array
191      */
192     private static PolynomialFunction buildPolynomial(final int degree,
193                                                       final ArrayList<BigFraction> coefficients,
194                                                       final RecurrenceCoefficientsGenerator generator) {
195 
196         final int maxDegree = (int) Math.floor(Math.sqrt(2 * coefficients.size())) - 1;
197         synchronized (PolynomialsUtils.class) {
198             if (degree > maxDegree) {
199                 computeUpToDegree(degree, maxDegree, generator, coefficients);
200             }
201         }
202 
203         // coefficient  for polynomial 0 is  l [0]
204         // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
205         // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
206         // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
207         // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
208         // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
209         // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
210         // ...
211         final int start = degree * (degree + 1) / 2;
212 
213         final double[] a = new double[degree + 1];
214         for (int i = 0; i <= degree; ++i) {
215             a[i] = coefficients.get(start + i).doubleValue();
216         }
217 
218         // build the polynomial
219         return new PolynomialFunction(a);
220 
221     }
222     
223     /** Compute polynomial coefficients up to a given degree.
224      * @param degree maximal degree
225      * @param maxDegree current maximal degree
226      * @param generator recurrence coefficients generator
227      * @param coefficients list where the computed coefficients should be appended
228      */
229     private static void computeUpToDegree(final int degree, final int maxDegree,
230                                           final RecurrenceCoefficientsGenerator generator,
231                                           final ArrayList<BigFraction> coefficients) {
232 
233         int startK = (maxDegree - 1) * maxDegree / 2;
234         for (int k = maxDegree; k < degree; ++k) {
235 
236             // start indices of two previous polynomials Pk(X) and Pk-1(X)
237             int startKm1 = startK;
238             startK += k;
239 
240             // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
241             BigFraction[] ai = generator.generate(k);
242 
243             BigFraction ck     = coefficients.get(startK);
244             BigFraction ckm1   = coefficients.get(startKm1);
245 
246             // degree 0 coefficient
247             coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
248 
249             // degree 1 to degree k-1 coefficients
250             for (int i = 1; i < k; ++i) {
251                 final BigFraction ckPrev = ck;
252                 ck     = coefficients.get(startK + i);
253                 ckm1   = coefficients.get(startKm1 + i);
254                 coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
255             }
256 
257             // degree k coefficient
258             final BigFraction ckPrev = ck;
259             ck = coefficients.get(startK + k);
260             coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
261 
262             // degree k+1 coefficient
263             coefficients.add(ck.multiply(ai[1]));
264 
265         }
266 
267     }
268 
269     /** Interface for recurrence coefficients generation. */
270     private static interface RecurrenceCoefficientsGenerator {
271         /**
272          * Generate recurrence coefficients.
273          * @param k highest degree of the polynomials used in the recurrence
274          * @return an array of three coefficients such that
275          * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
276          */
277         BigFraction[] generate(int k);
278     }
279 
280 }