This is almost certainly a consequence of numpy sometimes using pairwise
summation and sometimes
not.
Let’s build a diagnostic array:
eps = (np.nextafter(1.0, 2)-1.0) / 21+eps+eps+eps# 1.0(1+eps)+(eps+eps)# 1.0000000000000002X = np.full((32, 32), eps)X[0, 0] = 1X.sum(0)[0]# 1.0X.sum(1)[0]# 1.000000000000003X[:, 0].sum()# 1.000000000000003
This strongly suggests that 1D arrays and contiguous axes use pairwise
summation while strided axes in a multidimensional array don’t.
Note that to see that effect the array has to be large enough, otherwise numpy
falls back to ordinary summation.
