float32 np.cumsum is a naive running sum: 8.8% error at 1e7 terms

np.sum uses pairwise summation, np.cumsum (and np.add.accumulate, pandas Series.cumsum) cannot: each prefix is one rounded running sum, so the error grows with n. NumPy 2.5.4 / pandas 3.1.1, relative error of the last prefix vs math.fsum: float32 array of 0.1, n=1e5/1e6/1e7 -> 1.4e-4 / 9.6e-3 / 8.8e-2, while np.sum on the same array stays at ~1e-7. Uniform(0,1) float32 data: 2.8e-6 / 9.3e-6 / 2.8e-5 (sum: ~1e-8). Float64 cumsum of 1e7 U(0,1): 1.6e-14. The cumsum-difference moving average (cs[k:]-cs[:-k])/k on the 0.1 array returns 0.125 everywhere near the end, because once the running total passes ~1e6 the float32 spacing is 0.0625 or larger. Fix: np.cumsum(x, dtype=np.float64) (exact to printed precision here) or compute windowed sums with convolve or a rolling method instead of differencing a long prefix sum.