float32 np.dot is 80x less accurate than np.sum(a*b) at 1e7
For float32 vectors, np.dot / a@b go to BLAS sdot (scipy-openblas 0.3.34, NumPy 2.5.4). sdot accumulates in a few float32 SIMD lanes, so the error grows roughly linearly with n, while np.sum uses pairwise summation. Relative error vs an fsum reference (uniform [0,1) inputs): n=1e4: dot 6.8e-8, sum(a*b) 3.0e-8; n=1e6: 4.9e-7 vs 7.7e-9; n=1e7: 5.3e-6 vs 6.4e-8. The same with 1 or 8 BLAS threads. A naive running sum (cumsum[-1]) reached 1.2e-2 at 1e7. A float32 matvec with rows of 1e5 had max relative error 2.7e-7. If you need float32 storage but an accurate long reduction, use np.dot(a.astype(np.float64), b.astype(np.float64)) (2e-15), or np.sum(a*b) when memory for a temporary is fine. Expect sdot-level error from float32 matmul with long inner dimensions.
Author correction after re-running with 5 seeds per n (NumPy 2.5.4, OpenBLAS 0.3.34 SkylakeX kernel, 1 thread, median relative error of np.dot vs fsum): 'grows roughly linearly with n' is wrong. The single-seed n=1e6 figure in the post (4.9e-7) was the worst of the 5 seeds; the median is 6.3e-8. Medians: n=1e4 3.4e-8, 1e5 4.3e-8, 1e6 6.3e-8, 2e6 1.4e-7, 3e6 3.9e-7, 4e6 6.5e-7, 6e6 2.5e-6, 1e7 5.5e-6, 3e7 3.8e-5. So it is flat (no worse than np.sum(a*b), 1-3e-8 throughout) up to about 1e6, then grows faster than n (about 40x from 2e6 to 1e7). That fits biased rounding once the float32 lane accumulators get large. The headline (about 100x worse at 1e7) and the advice stand, and the breakpoint will depend on the BLAS kernel and SIMD width.