float32 statistics: np.var along an outer axis returned 17139 for true variance 1
NumPy 2.5.4, 1e7 float32 values 1e4 + N(0,1). Over the flat array np.var is fine (0.9995021 vs 0.99950196 in float64), because NumPy uses two passes plus pairwise summation. The textbook mean(x*x) - mean(x)**2 in float32 gives exactly 0.0. After reshaping to (1e6, 10), X.mean(axis=0) gives 10130.75 instead of 10000.00 and X.var(axis=0) gives 17139.6 instead of 0.9989. That axis is strided, so the accumulation is plain sequential float32 with no pairwise summation, and the ulp near 1e10 is 1024. Passing dtype=np.float64 to mean/var fixes it (0.99891), and so does reducing a contiguous copy. cumsum is always sequential: 2e7 float32 copies of 0.1 stall at 2**21 = 2097152 (0.1 is below half an ulp there), while .sum() gives 2000000.2. Rule: reduce float32 data with dtype=np.float64 unless you have checked the axis and the size.
Reproduced on NumPy 2.5.4 / CPython 3.12 with default_rng(0) instead of your data. Flat np.var is 0.999008 (float64: 0.999007926). Reshaped to (1e6, 10), axis=0 gives mean 10130.75 and var ≈17139, while dtype=np.float64 or a contiguous copy gives 0.998–0.999. cumsum of 2e7 float32 0.1s stops at 2097152 while .sum() is 2000000.2. One difference: with my data, the textbook mean(x*x) - mean(x)**2 in float32 came out as -8.0, not 0.0. The cancellation result depends on the data, and it can even be a negative 'variance', which makes a better canary than 0.0 because a later sqrt turns it into nan.