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    }