I agree with the sentiment that roundoff errors are expected when doing sums, but in this particular case I admit I’m mystified.
x + (-x) is always exactlly == 0 in floating-point arithmetic (for finite numbers), and similarly I think you should always have x * (-y) == -(x * y) exactly in the default rounding mode (round-to-nearest/ties-to-even). The combination of these two rules seems like it should mean that [x -x] * [y, y] == 0 exactly.
I’m not sure how you could implement the dot product to get a nonzero result here, short of changing the rounding mode? It’s a puzzle.
Mind you, I’m not saying it’s a bug — the -3.6e-16 answer is well within the tolerance I would generically expect for floating-point calculations on numbers of magnitude ≈ 4 — just that I’m curious how the BLAS dot could be implemented differently for length-2 vectors here. (I’m guessing it’s some SIMD thing, but I’m still not sure how SIMD would change the rounding in this case.)