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