Add hermite polynomial coefficients

meson.build

@@ -46,6 +46,7 @@ sources_list = [
     'mir/interpolate/utility',
     'mir/lob',
     'mir/math/func/expdigamma',
+    'mir/math/func/hermite',
     'mir/math/func/normal',
     'mir/math/numeric',
     'mir/math/stat',

source/mir/math/func/hermite.d

@@ -0,0 +1,182 @@
+/++
+Hermite Polynomial Coefficients
+
+License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0)
+Authors: John Hall
+
++/
+
+module mir.math.func.hermite;
+
+/++
+Normalized (probabilist's) Hermite polynomial coefficients
+
+Params:
+    N = Degree of polynomial
+
+See_also:
+    $(LINK2 https://en.wikipedia.org/wiki/Hermite_polynomials, Hermite polynomials)
++/
+@safe pure nothrow
+long[] hermiteCoefficientsNorm(size_t N)
+{
+    if (N == 0) {
+      return [1];  
+    } else {
+        typeof(return) output = [0, 1];
+        if (N > 1) {
+            output.length = N + 1;
+            int K;
+            typeof(return) h_N_minus_1 = hermiteCoefficientsNorm(0); // to be copied to h_N_minus_2 in loop
+            h_N_minus_1.length = N;
+            typeof(return) h_N_minus_2; //value replaced later
+            h_N_minus_2.length = N - 1;
+
+            foreach (size_t j; 2..(N + 1)) {
+                h_N_minus_2[0..(j - 1)] = h_N_minus_1[0..(j - 1)];
+                h_N_minus_1[0..j] = output[0..j];
+                K = -(cast(int) j - 1);
+                output[0] = K * h_N_minus_2[0];
+                foreach (size_t i; 1..(j - 1)) {
+                    output[i] = h_N_minus_1[i - 1] + K * h_N_minus_2[i];
+                }
+                foreach (size_t i; (j - 1)..(j + 1)) {
+                    output[i] = h_N_minus_1[i - 1];
+                }
+            }
+        }
+        return output;
+    }
+}
+
+///
+version(mir_test)
+@safe pure nothrow
+unittest
+{
+    import mir.polynomial: polynomial;
+    import mir.rc.array: rcarray;
+    import mir.test: should;
+
+    auto h0 = hermiteCoefficientsNorm(0).rcarray!(const double).polynomial;
+    auto h1 = hermiteCoefficientsNorm(1).rcarray!(const double).polynomial;
+    auto h2 = hermiteCoefficientsNorm(2).rcarray!(const double).polynomial;
+    auto h3 = hermiteCoefficientsNorm(3).rcarray!(const double).polynomial;
+    auto h4 = hermiteCoefficientsNorm(4).rcarray!(const double).polynomial;
+    auto h5 = hermiteCoefficientsNorm(5).rcarray!(const double).polynomial;
+    auto h6 = hermiteCoefficientsNorm(6).rcarray!(const double).polynomial;
+    auto h7 = hermiteCoefficientsNorm(7).rcarray!(const double).polynomial;
+    auto h8 = hermiteCoefficientsNorm(8).rcarray!(const double).polynomial;
+    auto h9 = hermiteCoefficientsNorm(9).rcarray!(const double).polynomial;
+    auto h10 = hermiteCoefficientsNorm(10).rcarray!(const double).polynomial;
+
+    h0(3).should == 1;
+    h1(3).should == 3;
+    h2(3).should == 8;
+    h3(3).should == 18;
+    h4(3).should == 30;
+    h5(3).should == 18;
+    h6(3).should == -96;
+    h7(3).should == -396;
+    h8(3).should == -516;
+    h9(3).should == 1_620;
+    h10(3).should == 9_504;
+}
+
+/// Also works with @nogc CTFE
+version(mir_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import mir.ndslice.slice: sliced;
+    import mir.test: should;
+
+    static immutable result = [-1, 0, 1];
+
+    static immutable hc2 = hermiteCoefficientsNorm(2);
+    hc2.sliced.should == result;
+}
+
+/++
+Physicist's Hermite polynomial coefficients
+
+Params:
+    N = Degree of polynomial
+
+See_also:
+    $(LINK2 https://en.wikipedia.org/wiki/Hermite_polynomials, Hermite polynomials)
++/
+@safe pure nothrow
+long[] hermiteCoefficients(size_t N)
+{
+    if (N == 0) {
+      return [1];  
+    } else {
+        typeof(return) output = [0, 2];
+        if (N > 1) {
+            output.length = N + 1;
+            typeof(return) h_N_minus_1 = hermiteCoefficients(0);
+            h_N_minus_1.length = N;
+
+            foreach (size_t j; 2..(N + 1)) {
+                h_N_minus_1[0..j] = output[0..j];
+                output[0] = -h_N_minus_1[1];
+                foreach (size_t i; 1..(j - 1)) {
+                    output[i] = 2 * h_N_minus_1[i - 1] - (cast(int) i + 1) * h_N_minus_1[i + 1];
+                }
+                foreach (size_t i; (j - 1)..(j + 1)) {
+                    output[i] = 2 * h_N_minus_1[i - 1];
+                }
+            }
+        }
+        return output;
+    }
+}
+
+///
+version(mir_test)
+@safe pure nothrow
+unittest
+{
+    import mir.polynomial: polynomial;
+    import mir.rc.array: rcarray;
+    import mir.test: should;
+
+    auto h0 = hermiteCoefficients(0).rcarray!(const double).polynomial;
+    auto h1 = hermiteCoefficients(1).rcarray!(const double).polynomial;
+    auto h2 = hermiteCoefficients(2).rcarray!(const double).polynomial;
+    auto h3 = hermiteCoefficients(3).rcarray!(const double).polynomial;
+    auto h4 = hermiteCoefficients(4).rcarray!(const double).polynomial;
+    auto h5 = hermiteCoefficients(5).rcarray!(const double).polynomial;
+    auto h6 = hermiteCoefficients(6).rcarray!(const double).polynomial;
+    auto h7 = hermiteCoefficients(7).rcarray!(const double).polynomial;
+    auto h8 = hermiteCoefficients(8).rcarray!(const double).polynomial;
+    auto h9 = hermiteCoefficients(9).rcarray!(const double).polynomial;
+    auto h10 = hermiteCoefficients(10).rcarray!(const double).polynomial;
+
+    h0(3).should == 1;
+    h1(3).should == 6;
+    h2(3).should == 34;
+    h3(3).should == 180;
+    h4(3).should == 876;
+    h5(3).should == 3_816;
+    h6(3).should == 14_136;
+    h7(3).should == 39_024;
+    h8(3).should == 36_240;
+    h9(3).should == -406_944;
+    h10(3).should == -3_093_984;
+}
+
+/// Also works with @nogc CTFE
+version(mir_test)
+@safe pure nothrow @nogc
+unittest
+{
+    import mir.ndslice.slice: sliced;
+    import mir.test: should;
+
+    static immutable result = [-2, 0, 4];
+
+    static immutable hc2 = hermiteCoefficients(2);
+    hc2.sliced.should == result;
+}