np.polymul return type is np.float64 or np.complex128 when given an all-zero argument

[roryyorke]

np.polymul returns an object with type np.float64 when one argument is all zero, and both arguments have type np.int64 or np.float32. Something similar happens with all zero np.complex64 giving result type np.complex128.

This doesn't happen with non-zero arguments; there the result is as expected.

This bug isn't present in np.convolve.

Reproducing code example:

>>> import numpy as np
>>> np.__version__
'1.18.4'
>>> a = np.array([1,2,3])
>>> z = np.array([0,0,0])
>>> np.polymul(a.astype(np.int64), a.astype(np.int64)).dtype
dtype('int64')
>>> np.polymul(a.astype(np.int64), z.astype(np.int64)).dtype
dtype('float64')
>>> np.polymul(a.astype(np.float32), z.astype(np.float32)).dtype
dtype('float64')
>>> np.polymul(a.astype(np.complex64), z.astype(np.complex64)).dtype
dtype('complex128')

Numpy/Python version information:

>>> import sys, numpy; print(numpy.__version__, sys.version)
1.18.4 3.7.5 (default, Nov  7 2019, 10:50:52) 
[GCC 8.3.0]

[rossbar]

This happens because polymul converts the input arrays to poly1d instances under-the-hood. The poly1d constructor has a special catch for all-zero coefficient arrays:

numpy/numpy/lib/polynomial.py

Lines 1170 to 1172 in ba5205c

	c_or_r = trim_zeros(c_or_r, trim='f')
	if len(c_or_r) == 0:
	c_or_r = NX.array([0.])

All that said, it's not clear what should be done about it. It would be straightforward to make np.array([0]) (see lines above) conform to the dtype of the input array, but this would be a backwards-incompatible change. Generally, user's are encouraged to use the newer np.polynomial, particularly the class-based interface (e.g. np.polynomial.Polynomial) which should (in principle) have more consistent behavior.

[eric-wieser]

I'd say this one is a clear bug that we should just fix - but I agree with the rest of your comment.