Function hammersley_sphere 

pub fn hammersley_sphere(n: usize, k: usize) -> Vec<Vec<f64>>

Hammersley point sequence mapped to the unit (n − 1)-sphere via inverse-normal-CDF + normalize (Marsaglia method on a low-discrepancy sequence). Returns k row vectors of length n, each unit-norm and asymptotically uniform on the sphere for any n as k → ∞.

Construction

n ≥ 3 (Marsaglia): for each i = 1..=k we generate n coordinates in (0, 1):

• u_0[i] = (i − 0.5) / k (linear, half-shifted to avoid endpoint bias),
• u_d[i] = radical_inverse(i, primes[d − 1]) for d = 1..n (van der Corput in the first n − 1 primes).

Each coordinate is mapped to a standard-normal sample via inverse CDF (Φ⁻¹), then the resulting n-vector is normalized to unit length. By the rotational invariance of the multivariate normal distribution, this produces points uniformly distributed on the unit sphere for any n, and the smooth Φ⁻¹ mapping makes the resulting sphere point set quasi-random — well-distributed in practice, though formally “low-discrepancy on the sphere” requires bounds (Brauchart & Dick on QMC on spheres) beyond Koksma-Hlawka because Φ⁻¹ is non-Lipschitz near 0 and 1.

n = 2 (half-circle): a Marsaglia draw on a 2-sphere would cover the full unit circle, but in MEMD a projection onto v and onto −v produces sign-flipped extrema sets (max ↔ min) and an identical per-channel mean-envelope contribution. The full-circle set would therefore contain antipodal redundancy at 2× the necessary compute. For n = 2 we instead lay k directions uniformly on the upper semicircle (θ_i = π·(i − 0.5)/k, v_i = (cos θ_i, sin θ_i)), preserving the antipode-redundancy optimization.

Why not direct hyperspherical coordinates for n ≥ 3?

The naive θ_1 = arccos(1 − 2u_1), θ_d = 2π·u_d mapping is the correct marginal-CDF inverse only for n = 3: the joint density on hyperspherical angles is ∝ sin^{n−2}(θ_1) · sin^{n−3}(θ_2) · … · sin(θ_{n−2}), which lacks closed-form CDF inverses for n ≥ 4. Earlier versions of this code used the naive form and produced biased direction sets for higher-dimensional MEMD inputs. The Marsaglia/Gaussian-normalize approach is the standard robust alternative.

References

• Marsaglia, G. (1972). Choosing a Point from the Surface of a Sphere. Ann. Math. Statist. 43(2), 645–646.
• Fang, K.-T., & Wang, Y. (1994). Number-Theoretic Methods in Statistics. Chapman & Hall — Hammersley/Halton point construction.
• Brauchart, J. S., & Dick, J. (2013). Quasi-Monte Carlo rules for numerical integration over the unit sphere. Numer. Math. 121(3), 473–502 — on sphere QMC discrepancy beyond Koksma-Hlawka.