np.vdot() incorrect on signed zeros

[pitrou]

See:

>>> complex(1,0) * complex(-1,0).conjugate()
(-1-0j)

but:

>>> a = np.array([complex(1, 0)])
>>> b = np.array([complex(-1, 0)])
>>> np.vdot(a, b)
(-1+0j)

Also:

>>> a * b.conjugate()
array([-1.-0.j])

[pv]

I'd be surprised if BLAS implementations are consistent vs. signed zero, there's probably a speed vs ieee-correctness tradeoff. vdot appears to dispatch to cblas, so it may be the inconsistency comes from there.

[pitrou]

For contiguous arrays, it probably depends on the BLAS implementation. But the non-contiguous implementation is faulty:

>>> a = np.array([complex(1, 0)]*2)[::-1]
>>> b = np.array([complex(-1, 0)]*2)[::-1]
>>> (a * b.conjugate())
array([-1.-0.j, -1.-0.j])
>>> np.vdot(a, b)
(-2+0j)

Note that np.sum() is faulty as well:

>>> (a * b.conjugate()).sum()
(-2+0j)

[pitrou]

By the way, the way this popped up is that MKL (which is becoming the default BLAS in Anaconda) seems to get it right.

[pitrou]

Though it may depend on the particular specialization (depending on the ISA features)... Perhaps we should skip testing for negative zeroes ;)

[seberg]

Pretty much the same as gh-21342, so will just make a note there. IIRC from the PR: You can possibly fix this for multi-dimenasional dot (0 initial value is passed in and -0 might help), but not vdot, as the scalar version doesn't pass it in: Probably not worthwhile, one issue is enough.