Function gauss_hermite_impl 

pub fn gauss_hermite_impl<P, E, V>(k: V, n: u32) -> (V, V)
where
    P: Policy,
    E: FloatElement,
    V: FloatVector<Element = E> + SpecializedSpecialMath<E>,

The k-th node and weight of the n-point Gauss-Hermite rule (weight e^{-x^2} on the line), k per lane, k = 0 the largest root.

Seed: the WKB phase of the Hermite equation, x = sqrt(2n+1) cos(phi) with phi - sin(2 phi)/2 = 2 pi (k + 3/4)/(2n + 1), solved per lane by four Newton steps from phi = (3c/2)^{1/3}. At the edge this reproduces the Airy constant (x_0 ~ sqrt(2n+1) - 1.856 (2n+1)^{-1/6}) to three digits. The left half is the mirror of the right. Newton then runs on h_m = H_m / m!, whose recurrence h_{m+1} = (2x h_m - 2 h_{m-1})/(m + 1) has a scalar divisor and stays in range where the raw H_m overflows at degree 48, with h_n' = 2 h_{n-1}. Weight w = sqrt(pi) 2^{n-1} / ((n-1)! n h_{n-1}(x_k)^2), the scalar factor a running product. That factor underflows past n = 170 in f64 and n = 40 in f32, which bounds the rule.