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    
023    /**
024     * A collection of static methods that operate on or return polynomials.
025     *
026     * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
027     * @since 2.0
028     */
029    public class PolynomialsUtils {
030    
031        /** Coefficients for Chebyshev polynomials. */
032        private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;
033    
034        /** Coefficients for Hermite polynomials. */
035        private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;
036    
037        /** Coefficients for Laguerre polynomials. */
038        private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;
039    
040        /** Coefficients for Legendre polynomials. */
041        private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;
042    
043        static {
044    
045            // initialize recurrence for Chebyshev polynomials
046            // T0(X) = 1, T1(X) = 0 + 1 * X
047            CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
048            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
049            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
050            CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
051    
052            // initialize recurrence for Hermite polynomials
053            // H0(X) = 1, H1(X) = 0 + 2 * X
054            HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
055            HERMITE_COEFFICIENTS.add(BigFraction.ONE);
056            HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
057            HERMITE_COEFFICIENTS.add(BigFraction.TWO);
058    
059            // initialize recurrence for Laguerre polynomials
060            // L0(X) = 1, L1(X) = 1 - 1 * X
061            LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
062            LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
063            LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
064            LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
065    
066            // initialize recurrence for Legendre polynomials
067            // P0(X) = 1, P1(X) = 0 + 1 * X
068            LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
069            LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
070            LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
071            LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
072    
073        }
074    
075        /**
076         * Private constructor, to prevent instantiation.
077         */
078        private PolynomialsUtils() {
079        }
080    
081        /**
082         * Create a Chebyshev polynomial of the first kind.
083         * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
084         * polynomials of the first kind</a> are orthogonal polynomials.
085         * They can be defined by the following recurrence relations:
086         * <pre>
087         *  T<sub>0</sub>(X)   = 1
088         *  T<sub>1</sub>(X)   = X
089         *  T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
090         * </pre></p>
091         * @param degree degree of the polynomial
092         * @return Chebyshev polynomial of specified degree
093         */
094        public static PolynomialFunction createChebyshevPolynomial(final int degree) {
095            return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
096                    new RecurrenceCoefficientsGenerator() {
097                private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
098                /** {@inheritDoc} */
099                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    }