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 */
017package org.apache.commons.math3.analysis.polynomials;
018
019import java.util.ArrayList;
020import java.util.HashMap;
021import java.util.List;
022import java.util.Map;
023
024import org.apache.commons.math3.fraction.BigFraction;
025import org.apache.commons.math3.util.CombinatoricsUtils;
026import org.apache.commons.math3.util.FastMath;
027
028/**
029 * A collection of static methods that operate on or return polynomials.
030 *
031 * @version $Id: PolynomialsUtils.java 1517203 2013-08-24 21:55:35Z psteitz $
032 * @since 2.0
033 */
034public class PolynomialsUtils {
035
036    /** Coefficients for Chebyshev polynomials. */
037    private static final List<BigFraction> CHEBYSHEV_COEFFICIENTS;
038
039    /** Coefficients for Hermite polynomials. */
040    private static final List<BigFraction> HERMITE_COEFFICIENTS;
041
042    /** Coefficients for Laguerre polynomials. */
043    private static final List<BigFraction> LAGUERRE_COEFFICIENTS;
044
045    /** Coefficients for Legendre polynomials. */
046    private static final List<BigFraction> LEGENDRE_COEFFICIENTS;
047
048    /** Coefficients for Jacobi polynomials. */
049    private static final Map<JacobiKey, List<BigFraction>> JACOBI_COEFFICIENTS;
050
051    static {
052
053        // initialize recurrence for Chebyshev polynomials
054        // T0(X) = 1, T1(X) = 0 + 1 * X
055        CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
056        CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
057        CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
058        CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
059
060        // initialize recurrence for Hermite polynomials
061        // H0(X) = 1, H1(X) = 0 + 2 * X
062        HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
063        HERMITE_COEFFICIENTS.add(BigFraction.ONE);
064        HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
065        HERMITE_COEFFICIENTS.add(BigFraction.TWO);
066
067        // initialize recurrence for Laguerre polynomials
068        // L0(X) = 1, L1(X) = 1 - 1 * X
069        LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
070        LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
071        LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
072        LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);
073
074        // initialize recurrence for Legendre polynomials
075        // P0(X) = 1, P1(X) = 0 + 1 * X
076        LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
077        LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
078        LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
079        LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
080
081        // initialize map for Jacobi polynomials
082        JACOBI_COEFFICIENTS = new HashMap<JacobiKey, List<BigFraction>>();
083
084    }
085
086    /**
087     * Private constructor, to prevent instantiation.
088     */
089    private PolynomialsUtils() {
090    }
091
092    /**
093     * Create a Chebyshev polynomial of the first kind.
094     * <p><a href="http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
095     * polynomials of the first kind</a> are orthogonal polynomials.
096     * They can be defined by the following recurrence relations:
097     * <pre>
098     *  T<sub>0</sub>(X)   = 1
099     *  T<sub>1</sub>(X)   = X
100     *  T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
101     * </pre></p>
102     * @param degree degree of the polynomial
103     * @return Chebyshev polynomial of specified degree
104     */
105    public static PolynomialFunction createChebyshevPolynomial(final int degree) {
106        return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS,
107                new RecurrenceCoefficientsGenerator() {
108            private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
109            /** {@inheritDoc} */
110            public BigFraction[] generate(int k) {
111                return coeffs;
112            }
113        });
114    }
115
116    /**
117     * Create a Hermite polynomial.
118     * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
119     * polynomials</a> are orthogonal polynomials.
120     * They can be defined by the following recurrence relations:
121     * <pre>
122     *  H<sub>0</sub>(X)   = 1
123     *  H<sub>1</sub>(X)   = 2X
124     *  H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
125     * </pre></p>
126
127     * @param degree degree of the polynomial
128     * @return Hermite polynomial of specified degree
129     */
130    public static PolynomialFunction createHermitePolynomial(final int degree) {
131        return buildPolynomial(degree, HERMITE_COEFFICIENTS,
132                new RecurrenceCoefficientsGenerator() {
133            /** {@inheritDoc} */
134            public BigFraction[] generate(int k) {
135                return new BigFraction[] {
136                        BigFraction.ZERO,
137                        BigFraction.TWO,
138                        new BigFraction(2 * k)};
139            }
140        });
141    }
142
143    /**
144     * Create a Laguerre polynomial.
145     * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
146     * polynomials</a> are orthogonal polynomials.
147     * They can be defined by the following recurrence relations:
148     * <pre>
149     *        L<sub>0</sub>(X)   = 1
150     *        L<sub>1</sub>(X)   = 1 - X
151     *  (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
152     * </pre></p>
153     * @param degree degree of the polynomial
154     * @return Laguerre polynomial of specified degree
155     */
156    public static PolynomialFunction createLaguerrePolynomial(final int degree) {
157        return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
158                new RecurrenceCoefficientsGenerator() {
159            /** {@inheritDoc} */
160            public BigFraction[] generate(int k) {
161                final int kP1 = k + 1;
162                return new BigFraction[] {
163                        new BigFraction(2 * k + 1, kP1),
164                        new BigFraction(-1, kP1),
165                        new BigFraction(k, kP1)};
166            }
167        });
168    }
169
170    /**
171     * Create a Legendre polynomial.
172     * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
173     * polynomials</a> are orthogonal polynomials.
174     * They can be defined by the following recurrence relations:
175     * <pre>
176     *        P<sub>0</sub>(X)   = 1
177     *        P<sub>1</sub>(X)   = X
178     *  (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
179     * </pre></p>
180     * @param degree degree of the polynomial
181     * @return Legendre polynomial of specified degree
182     */
183    public static PolynomialFunction createLegendrePolynomial(final int degree) {
184        return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
185                               new RecurrenceCoefficientsGenerator() {
186            /** {@inheritDoc} */
187            public BigFraction[] generate(int k) {
188                final int kP1 = k + 1;
189                return new BigFraction[] {
190                        BigFraction.ZERO,
191                        new BigFraction(k + kP1, kP1),
192                        new BigFraction(k, kP1)};
193            }
194        });
195    }
196
197    /**
198     * Create a Jacobi polynomial.
199     * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi
200     * polynomials</a> are orthogonal polynomials.
201     * They can be defined by the following recurrence relations:
202     * <pre>
203     *        P<sub>0</sub><sup>vw</sup>(X)   = 1
204     *        P<sub>-1</sub><sup>vw</sup>(X)  = 0
205     *  2k(k + v + w)(2k + v + w - 2) P<sub>k</sub><sup>vw</sup>(X) =
206     *  (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) X + v<sup>2</sup> - w<sup>2</sup>] P<sub>k-1</sub><sup>vw</sup>(X)
207     *  - 2(k + v - 1)(k + w - 1)(2k + v + w) P<sub>k-2</sub><sup>vw</sup>(X)
208     * </pre></p>
209     * @param degree degree of the polynomial
210     * @param v first exponent
211     * @param w second exponent
212     * @return Jacobi polynomial of specified degree
213     */
214    public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) {
215
216        // select the appropriate list
217        final JacobiKey key = new JacobiKey(v, w);
218
219        if (!JACOBI_COEFFICIENTS.containsKey(key)) {
220
221            // allocate a new list for v, w
222            final List<BigFraction> list = new ArrayList<BigFraction>();
223            JACOBI_COEFFICIENTS.put(key, list);
224
225            // Pv,w,0(x) = 1;
226            list.add(BigFraction.ONE);
227
228            // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2
229            list.add(new BigFraction(v - w, 2));
230            list.add(new BigFraction(2 + v + w, 2));
231
232        }
233
234        return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key),
235                               new RecurrenceCoefficientsGenerator() {
236            /** {@inheritDoc} */
237            public BigFraction[] generate(int k) {
238                k++;
239                final int kvw      = k + v + w;
240                final int twoKvw   = kvw + k;
241                final int twoKvwM1 = twoKvw - 1;
242                final int twoKvwM2 = twoKvw - 2;
243                final int den      = 2 * k *  kvw * twoKvwM2;
244
245                return new BigFraction[] {
246                    new BigFraction(twoKvwM1 * (v * v - w * w), den),
247                    new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den),
248                    new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den)
249                };
250            }
251        });
252
253    }
254
255    /** Inner class for Jacobi polynomials keys. */
256    private static class JacobiKey {
257
258        /** First exponent. */
259        private final int v;
260
261        /** Second exponent. */
262        private final int w;
263
264        /** Simple constructor.
265         * @param v first exponent
266         * @param w second exponent
267         */
268        public JacobiKey(final int v, final int w) {
269            this.v = v;
270            this.w = w;
271        }
272
273        /** Get hash code.
274         * @return hash code
275         */
276        @Override
277        public int hashCode() {
278            return (v << 16) ^ w;
279        }
280
281        /** Check if the instance represent the same key as another instance.
282         * @param key other key
283         * @return true if the instance and the other key refer to the same polynomial
284         */
285        @Override
286        public boolean equals(final Object key) {
287
288            if ((key == null) || !(key instanceof JacobiKey)) {
289                return false;
290            }
291
292            final JacobiKey otherK = (JacobiKey) key;
293            return (v == otherK.v) && (w == otherK.w);
294
295        }
296    }
297
298    /**
299     * Compute the coefficients of the polynomial <code>P<sub>s</sub>(x)</code>
300     * whose values at point {@code x} will be the same as the those from the
301     * original polynomial <code>P(x)</code> when computed at {@code x + shift}.
302     * Thus, if <code>P(x) = &Sigma;<sub>i</sub> a<sub>i</sub> x<sup>i</sup></code>,
303     * then
304     * <pre>
305     *  <table>
306     *   <tr>
307     *    <td><code>P<sub>s</sub>(x)</td>
308     *    <td>= &Sigma;<sub>i</sub> b<sub>i</sub> x<sup>i</sup></code></td>
309     *   </tr>
310     *   <tr>
311     *    <td></td>
312     *    <td>= &Sigma;<sub>i</sub> a<sub>i</sub> (x + shift)<sup>i</sup></code></td>
313     *   </tr>
314     *  </table>
315     * </pre>
316     *
317     * @param coefficients Coefficients of the original polynomial.
318     * @param shift Shift value.
319     * @return the coefficients <code>b<sub>i</sub></code> of the shifted
320     * polynomial.
321     */
322    public static double[] shift(final double[] coefficients,
323                                 final double shift) {
324        final int dp1 = coefficients.length;
325        final double[] newCoefficients = new double[dp1];
326
327        // Pascal triangle.
328        final int[][] coeff = new int[dp1][dp1];
329        for (int i = 0; i < dp1; i++){
330            for(int j = 0; j <= i; j++){
331                coeff[i][j] = (int) CombinatoricsUtils.binomialCoefficient(i, j);
332            }
333        }
334
335        // First polynomial coefficient.
336        for (int i = 0; i < dp1; i++){
337            newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i);
338        }
339
340        // Superior order.
341        final int d = dp1 - 1;
342        for (int i = 0; i < d; i++) {
343            for (int j = i; j < d; j++){
344                newCoefficients[i + 1] += coeff[j + 1][j - i] *
345                    coefficients[j + 1] * FastMath.pow(shift, j - i);
346            }
347        }
348
349        return newCoefficients;
350    }
351
352
353    /** Get the coefficients array for a given degree.
354     * @param degree degree of the polynomial
355     * @param coefficients list where the computed coefficients are stored
356     * @param generator recurrence coefficients generator
357     * @return coefficients array
358     */
359    private static PolynomialFunction buildPolynomial(final int degree,
360                                                      final List<BigFraction> coefficients,
361                                                      final RecurrenceCoefficientsGenerator generator) {
362
363        final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
364        synchronized (PolynomialsUtils.class) {
365            if (degree > maxDegree) {
366                computeUpToDegree(degree, maxDegree, generator, coefficients);
367            }
368        }
369
370        // coefficient  for polynomial 0 is  l [0]
371        // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
372        // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
373        // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
374        // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
375        // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
376        // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
377        // ...
378        final int start = degree * (degree + 1) / 2;
379
380        final double[] a = new double[degree + 1];
381        for (int i = 0; i <= degree; ++i) {
382            a[i] = coefficients.get(start + i).doubleValue();
383        }
384
385        // build the polynomial
386        return new PolynomialFunction(a);
387
388    }
389
390    /** Compute polynomial coefficients up to a given degree.
391     * @param degree maximal degree
392     * @param maxDegree current maximal degree
393     * @param generator recurrence coefficients generator
394     * @param coefficients list where the computed coefficients should be appended
395     */
396    private static void computeUpToDegree(final int degree, final int maxDegree,
397                                          final RecurrenceCoefficientsGenerator generator,
398                                          final List<BigFraction> coefficients) {
399
400        int startK = (maxDegree - 1) * maxDegree / 2;
401        for (int k = maxDegree; k < degree; ++k) {
402
403            // start indices of two previous polynomials Pk(X) and Pk-1(X)
404            int startKm1 = startK;
405            startK += k;
406
407            // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
408            BigFraction[] ai = generator.generate(k);
409
410            BigFraction ck     = coefficients.get(startK);
411            BigFraction ckm1   = coefficients.get(startKm1);
412
413            // degree 0 coefficient
414            coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
415
416            // degree 1 to degree k-1 coefficients
417            for (int i = 1; i < k; ++i) {
418                final BigFraction ckPrev = ck;
419                ck     = coefficients.get(startK + i);
420                ckm1   = coefficients.get(startKm1 + i);
421                coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
422            }
423
424            // degree k coefficient
425            final BigFraction ckPrev = ck;
426            ck = coefficients.get(startK + k);
427            coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
428
429            // degree k+1 coefficient
430            coefficients.add(ck.multiply(ai[1]));
431
432        }
433
434    }
435
436    /** Interface for recurrence coefficients generation. */
437    private interface RecurrenceCoefficientsGenerator {
438        /**
439         * Generate recurrence coefficients.
440         * @param k highest degree of the polynomials used in the recurrence
441         * @return an array of three coefficients such that
442         * P<sub>k+1</sub>(X) = (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
443         */
444        BigFraction[] generate(int k);
445    }
446
447}