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