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="https://en.wikipedia.org/wiki/Chebyshev_polynomials">Chebyshev
094     * polynomials of the first kind</a> are orthogonal polynomials.
095     * They can be defined by the following recurrence relations:</p><p>
096     * \(
097     *    T_0(x) = 1 \\
098     *    T_1(x) = x \\
099     *    T_{k+1}(x) = 2x T_k(x) - T_{k-1}(x)
100     * \)
101     * </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            /** Fixed recurrence coefficients. */
109            private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };
110            /** {@inheritDoc} */
111            public BigFraction[] generate(int k) {
112                return coeffs;
113            }
114        });
115    }
116
117    /**
118     * Create a Hermite polynomial.
119     * <p><a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
120     * polynomials</a> are orthogonal polynomials.
121     * They can be defined by the following recurrence relations:</p><p>
122     * \(
123     *  H_0(x) = 1 \\
124     *  H_1(x) = 2x \\
125     *  H_{k+1}(x) = 2x H_k(X) - 2k H_{k-1}(x)
126     * \)
127     * </p>
128
129     * @param degree degree of the polynomial
130     * @return Hermite polynomial of specified degree
131     */
132    public static PolynomialFunction createHermitePolynomial(final int degree) {
133        return buildPolynomial(degree, HERMITE_COEFFICIENTS,
134                new RecurrenceCoefficientsGenerator() {
135            /** {@inheritDoc} */
136            public BigFraction[] generate(int k) {
137                return new BigFraction[] {
138                        BigFraction.ZERO,
139                        BigFraction.TWO,
140                        new BigFraction(2 * k)};
141            }
142        });
143    }
144
145    /**
146     * Create a Laguerre polynomial.
147     * <p><a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
148     * polynomials</a> are orthogonal polynomials.
149     * They can be defined by the following recurrence relations:</p><p>
150     * \(
151     *   L_0(x) = 1 \\
152     *   L_1(x) = 1 - x \\
153     *   (k+1) L_{k+1}(x) = (2k + 1 - x) L_k(x) - k L_{k-1}(x)
154     * \)
155     * </p>
156     * @param degree degree of the polynomial
157     * @return Laguerre polynomial of specified degree
158     */
159    public static PolynomialFunction createLaguerrePolynomial(final int degree) {
160        return buildPolynomial(degree, LAGUERRE_COEFFICIENTS,
161                new RecurrenceCoefficientsGenerator() {
162            /** {@inheritDoc} */
163            public BigFraction[] generate(int k) {
164                final int kP1 = k + 1;
165                return new BigFraction[] {
166                        new BigFraction(2 * k + 1, kP1),
167                        new BigFraction(-1, kP1),
168                        new BigFraction(k, kP1)};
169            }
170        });
171    }
172
173    /**
174     * Create a Legendre polynomial.
175     * <p><a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
176     * polynomials</a> are orthogonal polynomials.
177     * They can be defined by the following recurrence relations:</p><p>
178     * \(
179     *   P_0(x) = 1 \\
180     *   P_1(x) = x \\
181     *   (k+1) P_{k+1}(x) = (2k+1) x P_k(x) - k P_{k-1}(x)
182     * \)
183     * </p>
184     * @param degree degree of the polynomial
185     * @return Legendre polynomial of specified degree
186     */
187    public static PolynomialFunction createLegendrePolynomial(final int degree) {
188        return buildPolynomial(degree, LEGENDRE_COEFFICIENTS,
189                               new RecurrenceCoefficientsGenerator() {
190            /** {@inheritDoc} */
191            public BigFraction[] generate(int k) {
192                final int kP1 = k + 1;
193                return new BigFraction[] {
194                        BigFraction.ZERO,
195                        new BigFraction(k + kP1, kP1),
196                        new BigFraction(k, kP1)};
197            }
198        });
199    }
200
201    /**
202     * Create a Jacobi polynomial.
203     * <p><a href="http://mathworld.wolfram.com/JacobiPolynomial.html">Jacobi
204     * polynomials</a> are orthogonal polynomials.
205     * They can be defined by the following recurrence relations:</p><p>
206     * \(
207     *    P_0^{vw}(x) = 1 \\
208     *    P_{-1}^{vw}(x) = 0 \\
209     *    2k(k + v + w)(2k + v + w - 2) P_k^{vw}(x) = \\
210     *    (2k + v + w - 1)[(2k + v + w)(2k + v + w - 2) x + v^2 - w^2] P_{k-1}^{vw}(x) \\
211     *  - 2(k + v - 1)(k + w - 1)(2k + v + w) P_{k-2}^{vw}(x)
212     * \)
213     * </p>
214     * @param degree degree of the polynomial
215     * @param v first exponent
216     * @param w second exponent
217     * @return Jacobi polynomial of specified degree
218     */
219    public static PolynomialFunction createJacobiPolynomial(final int degree, final int v, final int w) {
220
221        // select the appropriate list
222        final JacobiKey key = new JacobiKey(v, w);
223
224        if (!JACOBI_COEFFICIENTS.containsKey(key)) {
225
226            // allocate a new list for v, w
227            final List<BigFraction> list = new ArrayList<BigFraction>();
228            JACOBI_COEFFICIENTS.put(key, list);
229
230            // Pv,w,0(x) = 1;
231            list.add(BigFraction.ONE);
232
233            // P1(x) = (v - w) / 2 + (2 + v + w) * X / 2
234            list.add(new BigFraction(v - w, 2));
235            list.add(new BigFraction(2 + v + w, 2));
236
237        }
238
239        return buildPolynomial(degree, JACOBI_COEFFICIENTS.get(key),
240                               new RecurrenceCoefficientsGenerator() {
241            /** {@inheritDoc} */
242            public BigFraction[] generate(int k) {
243                k++;
244                final int kvw      = k + v + w;
245                final int twoKvw   = kvw + k;
246                final int twoKvwM1 = twoKvw - 1;
247                final int twoKvwM2 = twoKvw - 2;
248                final int den      = 2 * k *  kvw * twoKvwM2;
249
250                return new BigFraction[] {
251                    new BigFraction(twoKvwM1 * (v * v - w * w), den),
252                    new BigFraction(twoKvwM1 * twoKvw * twoKvwM2, den),
253                    new BigFraction(2 * (k + v - 1) * (k + w - 1) * twoKvw, den)
254                };
255            }
256        });
257
258    }
259
260    /** Inner class for Jacobi polynomials keys. */
261    private static class JacobiKey {
262
263        /** First exponent. */
264        private final int v;
265
266        /** Second exponent. */
267        private final int w;
268
269        /** Simple constructor.
270         * @param v first exponent
271         * @param w second exponent
272         */
273        JacobiKey(final int v, final int w) {
274            this.v = v;
275            this.w = w;
276        }
277
278        /** Get hash code.
279         * @return hash code
280         */
281        @Override
282        public int hashCode() {
283            return (v << 16) ^ w;
284        }
285
286        /** Check if the instance represent the same key as another instance.
287         * @param key other key
288         * @return true if the instance and the other key refer to the same polynomial
289         */
290        @Override
291        public boolean equals(final Object key) {
292
293            if ((key == null) || !(key instanceof JacobiKey)) {
294                return false;
295            }
296
297            final JacobiKey otherK = (JacobiKey) key;
298            return (v == otherK.v) && (w == otherK.w);
299
300        }
301    }
302
303    /**
304     * Compute the coefficients of the polynomial \(P_s(x)\)
305     * whose values at point {@code x} will be the same as the those from the
306     * original polynomial \(P(x)\) when computed at {@code x + shift}.
307     * <p>
308     * More precisely, let \(\Delta = \) {@code shift} and let
309     * \(P_s(x) = P(x + \Delta)\).  The returned array
310     * consists of the coefficients of \(P_s\).  So if \(a_0, ..., a_{n-1}\)
311     * are the coefficients of \(P\), then the returned array
312     * \(b_0, ..., b_{n-1}\) satisfies the identity
313     * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
314     *
315     * @param coefficients Coefficients of the original polynomial.
316     * @param shift Shift value.
317     * @return the coefficients \(b_i\) of the shifted
318     * polynomial.
319     */
320    public static double[] shift(final double[] coefficients,
321                                 final double shift) {
322        final int dp1 = coefficients.length;
323        final double[] newCoefficients = new double[dp1];
324
325        // Pascal triangle.
326        final int[][] coeff = new int[dp1][dp1];
327        for (int i = 0; i < dp1; i++){
328            for(int j = 0; j <= i; j++){
329                coeff[i][j] = (int) CombinatoricsUtils.binomialCoefficient(i, j);
330            }
331        }
332
333        // First polynomial coefficient.
334        for (int i = 0; i < dp1; i++){
335            newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i);
336        }
337
338        // Superior order.
339        final int d = dp1 - 1;
340        for (int i = 0; i < d; i++) {
341            for (int j = i; j < d; j++){
342                newCoefficients[i + 1] += coeff[j + 1][j - i] *
343                    coefficients[j + 1] * FastMath.pow(shift, j - i);
344            }
345        }
346
347        return newCoefficients;
348    }
349
350
351    /** Get the coefficients array for a given degree.
352     * @param degree degree of the polynomial
353     * @param coefficients list where the computed coefficients are stored
354     * @param generator recurrence coefficients generator
355     * @return coefficients array
356     */
357    private static PolynomialFunction buildPolynomial(final int degree,
358                                                      final List<BigFraction> coefficients,
359                                                      final RecurrenceCoefficientsGenerator generator) {
360        synchronized (coefficients) {
361            final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1;
362            if (degree > maxDegree) {
363                computeUpToDegree(degree, maxDegree, generator, coefficients);
364            }
365        }
366
367        // coefficient  for polynomial 0 is  l [0]
368        // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
369        // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
370        // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
371        // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
372        // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
373        // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
374        // ...
375        final int start = degree * (degree + 1) / 2;
376
377        final double[] a = new double[degree + 1];
378        for (int i = 0; i <= degree; ++i) {
379            a[i] = coefficients.get(start + i).doubleValue();
380        }
381
382        // build the polynomial
383        return new PolynomialFunction(a);
384
385    }
386
387    /** Compute polynomial coefficients up to a given degree.
388     * @param degree maximal degree
389     * @param maxDegree current maximal degree
390     * @param generator recurrence coefficients generator
391     * @param coefficients list where the computed coefficients should be appended
392     */
393    private static void computeUpToDegree(final int degree, final int maxDegree,
394                                          final RecurrenceCoefficientsGenerator generator,
395                                          final List<BigFraction> coefficients) {
396
397        int startK = (maxDegree - 1) * maxDegree / 2;
398        for (int k = maxDegree; k < degree; ++k) {
399
400            // start indices of two previous polynomials Pk(X) and Pk-1(X)
401            int startKm1 = startK;
402            startK += k;
403
404            // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
405            BigFraction[] ai = generator.generate(k);
406
407            BigFraction ck     = coefficients.get(startK);
408            BigFraction ckm1   = coefficients.get(startKm1);
409
410            // degree 0 coefficient
411            coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));
412
413            // degree 1 to degree k-1 coefficients
414            for (int i = 1; i < k; ++i) {
415                final BigFraction ckPrev = ck;
416                ck     = coefficients.get(startK + i);
417                ckm1   = coefficients.get(startKm1 + i);
418                coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
419            }
420
421            // degree k coefficient
422            final BigFraction ckPrev = ck;
423            ck = coefficients.get(startK + k);
424            coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));
425
426            // degree k+1 coefficient
427            coefficients.add(ck.multiply(ai[1]));
428
429        }
430
431    }
432
433    /** Interface for recurrence coefficients generation. */
434    private interface RecurrenceCoefficientsGenerator {
435        /**
436         * Generate recurrence coefficients.
437         * @param k highest degree of the polynomials used in the recurrence
438         * @return an array of three coefficients such that
439         * \( P_{k+1}(x) = (a[0] + a[1] x) P_k(x) - a[2] P_{k-1}(x) \)
440         */
441        BigFraction[] generate(int k);
442    }
443
444}