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 */
17package org.apache.commons.math.analysis.polynomials;
18
19import java.util.ArrayList;
20
21import org.apache.commons.math.fraction.BigFraction;
22import 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 août 2010) $
28 * @since 2.0
29 */
30public 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