np.sum is pairwise only along the memory-contiguous axis
NumPy uses pairwise summation only when the reduction runs along the fast (contiguous) axis. If you reduce over an outer axis, for example X.sum(axis=0) on a C-ordered array or axis=1 on a Fortran-ordered one, NumPy adds whole rows elementwise, which is plain recursive summation. Measured on 1e7 float32 uniforms, relative error of the total: contiguous 1D sum 1.2e-8; (1e6,10) array summed with axis=0: 5.9e-6; Fortran (2,5e6) array summed with axis=1: 2.1e-5; np.cumsum: 4.7e-6 (always sequential). A strided 1D view such as x[::2].sum() was still accurate (4e-8). For float32 data, either pass dtype=np.float64 (error 2e-16 here) or reduce along the contiguous axis. Don't assume np.sum has pairwise accuracy for every axis.
Confirmed on NumPy 2.5.4 with 1e7 float32 uniforms (seed 0), relative errors: contiguous 1D 3.8e-8; Fortran (2, 5e6) summed along axis=1, 1.7e-5; cumsum, 1.7e-5; x[::2].sum(), 3.2e-8; dtype=np.float64, 0. One measurement note for the (1e6, 10) axis=0 case: measure each column against its own fsum, not the total of the ten column sums. The per-column errors were mean 1.3e-5 and max 3.2e-5, while the error of their total was only 5.9e-7 because the column errors partly cancel. A contiguous copy (np.ascontiguousarray(X.T)[j].sum()) brings each column back to 5.5e-8. So the effect is, if anything, larger than a total-based figure suggests.