/* * Adapted from: * https://apache.googlesource.com/commons-math/+/d40d9b48aee58eea159eaa27ee14654fefd7af04/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialsUtils.java */ /* * Licensed to the Apache Software Foundation (ASF) under one or more * contributor license agreements. See the NOTICE file distributed with * this work for additional information regarding copyright ownership. * The ASF licenses this file to You under the Apache License, Version 2.0 * (the "License"); you may not use this file except in compliance with * the License. You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ import org.apache.commons.math3.analysis.polynomials.*; import org.apache.commons.math3.util.*; import java.util.ArrayList; import java.util.HashMap; import java.util.List; import java.util.Map; //import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; //import org.apache.commons.math3.analysis.polynomials.PolynomialsUtils; import org.apache.commons.math3.fraction.BigFraction; import org.apache.commons.math3.util.Combinatorics; import org.apache.commons.math3.util.FastMath; /** * A collection of static methods that operate on or return polynomials. * * @version $Id$ * @since 2.0 */ public class Hermite { /** Coefficients for Hermite polynomials. */ private static final List HERMITE_COEFFICIENTS; static { // initialize recurrence for Hermite polynomials // H0(X) = 1, H1(X) = 0 + 2 * X HERMITE_COEFFICIENTS = new ArrayList(); HERMITE_COEFFICIENTS.add(BigFraction.ONE); HERMITE_COEFFICIENTS.add(BigFraction.ZERO); HERMITE_COEFFICIENTS.add(BigFraction.TWO); } /** * Private constructor, to prevent instantiation. */ private Hermite() { } /** * Create a Hermite polynomial. *

Hermite * polynomials are orthogonal polynomials. * They can be defined by the following recurrence relations: *

         *  H0(X)   = 1
         *  H1(X)   = 2X
         *  Hk+1(X) = 2X Hk(X) - 2k Hk-1(X)
         * 

* @param degree degree of the polynomial * @return Hermite polynomial of specified degree */ public static PolynomialFunction createHermitePolynomial(final int degree) { return buildPolynomial(degree, HERMITE_COEFFICIENTS, new RecurrenceCoefficientsGenerator() { /** {@inheritDoc} */ public BigFraction[] generate(int k) { return new BigFraction[] { BigFraction.ZERO, BigFraction.TWO, new BigFraction(2 * k)}; } }); } /** * Compute the coefficients of the polynomial Ps(x) * whose values at point {@code x} will be the same as the those from the * original polynomial P(x) when computed at {@code x + shift}. * Thus, if P(x) = Σi ai xi, * then *
         *  
         *   
         *    
         *    
         *   
         *   
         *    
         *    
         *   
         *  
Ps(x)= Σi bi xi
= Σi ai (x + shift)i
*
* * @param coefficients Coefficients of the original polynomial. * @param shift Shift value. * @return the coefficients bi of the shifted * polynomial. */ public static double[] shift(final double[] coefficients, final double shift) { final int dp1 = coefficients.length; final double[] newCoefficients = new double[dp1]; // Pascal triangle. final int[][] coeff = new int[dp1][dp1]; for (int i = 0; i < dp1; i++){ for(int j = 0; j <= i; j++){ coeff[i][j] = (int) CombinatoricsUtils.binomialCoefficient(i, j); } } // First polynomial coefficient. for (int i = 0; i < dp1; i++){ newCoefficients[0] += coefficients[i] * FastMath.pow(shift, i); } // Superior order. final int d = dp1 - 1; for (int i = 0; i < d; i++) { for (int j = i; j < d; j++){ newCoefficients[i + 1] += coeff[j + 1][j - i] * coefficients[j + 1] * FastMath.pow(shift, j - i); } } return newCoefficients; } /** Get the coefficients array for a given degree. * @param degree degree of the polynomial * @param coefficients list where the computed coefficients are stored * @param generator recurrence coefficients generator * @return coefficients array */ private static PolynomialFunction buildPolynomial(final int degree, final List coefficients, final RecurrenceCoefficientsGenerator generator) { final int maxDegree = (int) FastMath.floor(FastMath.sqrt(2 * coefficients.size())) - 1; synchronized (PolynomialsUtils.class) { if (degree > maxDegree) { computeUpToDegree(degree, maxDegree, generator, coefficients); } } // coefficient for polynomial 0 is l [0] // coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1) // coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2) // coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3) // coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4) // coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5) // coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6) // ... final int start = degree * (degree + 1) / 2; final double[] a = new double[degree + 1]; for (int i = 0; i <= degree; ++i) { a[i] = coefficients.get(start + i).doubleValue(); } // build the polynomial return new PolynomialFunction(a); } /** Compute polynomial coefficients up to a given degree. * @param degree maximal degree * @param maxDegree current maximal degree * @param generator recurrence coefficients generator * @param coefficients list where the computed coefficients should be appended */ private static void computeUpToDegree(final int degree, final int maxDegree, final RecurrenceCoefficientsGenerator generator, final List coefficients) { int startK = (maxDegree - 1) * maxDegree / 2; for (int k = maxDegree; k < degree; ++k) { // start indices of two previous polynomials Pk(X) and Pk-1(X) int startKm1 = startK; startK += k; // Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X) BigFraction[] ai = generator.generate(k); BigFraction ck = coefficients.get(startK); BigFraction ckm1 = coefficients.get(startKm1); // degree 0 coefficient coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2]))); // degree 1 to degree k-1 coefficients for (int i = 1; i < k; ++i) { final BigFraction ckPrev = ck; ck = coefficients.get(startK + i); ckm1 = coefficients.get(startKm1 + i); coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2]))); } // degree k coefficient final BigFraction ckPrev = ck; ck = coefficients.get(startK + k); coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1]))); // degree k+1 coefficient coefficients.add(ck.multiply(ai[1])); } } /** Interface for recurrence coefficients generation. */ private interface RecurrenceCoefficientsGenerator { /** * Generate recurrence coefficients. * @param k highest degree of the polynomials used in the recurrence * @return an array of three coefficients such that * Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X) */ BigFraction[] generate(int k); } }