Repository: numpy/numpy
Pull request: BUG: Respect dtype of all-zero argument to poly1d
Fixes gh-16354. 
Previously, `np.poly1d(z).coeffs.dtype` would always be `np.float64` for zero array `z`, regardless of `z`'s dtype.

This led to unexpected dtype promotion in methods such as `np.polymul`:
```python
>>> np.polymul(np.array([1, 2, 3], dtype=np.int64), np.array([0, 0, 0], dtype=np.int64)).dtype
dtype('float64')
```

This PR retains the dtype of zero array inputs.
```python
>>> np.polymul(np.array([1, 2, 3], dtype=np.int64), np.array([0, 0, 0], dtype=np.int64)).dtype
dtype('int64')
```


Problem statement:
np.polymul return type is np.float64 or np.complex128 when given an all-zero argument
`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:
```python
>>> 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]
```

Hints from issue discussion:
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:

https://github.com/numpy/numpy/blob/ba5205c1b8a7bebdf37535c701ffcd8bd79362bb/numpy/lib/polynomial.py#L1170-L1172

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.
I'd say this one is a clear bug that we should just fix - but I agree with the rest of your comment.

Commit message: BUG: Respect dtype of all-zero argument to poly1d

Fixes gh-16354. Previously np.poly1d(z).coeffs.dtype would always
be np.float64 for zero array z, regardless of z's dtype.

Patch to review:
diff --git a/numpy/lib/polynomial.py b/numpy/lib/polynomial.py
index 7b89eeb70367..28c3f5e99cd8 100644
--- a/numpy/lib/polynomial.py
+++ b/numpy/lib/polynomial.py
@@ -1236,7 +1236,7 @@ def __init__(self, c_or_r, r=False, variable=None):
             raise ValueError("Polynomial must be 1d only.")
         c_or_r = trim_zeros(c_or_r, trim='f')
         if len(c_or_r) == 0:
-            c_or_r = NX.array([0.])
+            c_or_r = NX.array([0.], dtype=c_or_r.dtype)
         self._coeffs = c_or_r
         if variable is None:
             variable = 'x'
diff --git a/numpy/lib/tests/test_polynomial.py b/numpy/lib/tests/test_polynomial.py
index ab6691b4384f..6c3e4fa02212 100644
--- a/numpy/lib/tests/test_polynomial.py
+++ b/numpy/lib/tests/test_polynomial.py
@@ -227,6 +227,20 @@ def test_poly_int_overflow(self):
         v = np.arange(1, 21)
         assert_almost_equal(np.poly(v), np.poly(np.diag(v)))
 
+    def test_zero_poly_dtype(self):
+        """
+        Regression test for gh-16354.
+        """
+        z = np.array([0, 0, 0])
+        p = np.poly1d(z.astype(np.int64))
+        assert_equal(p.coeffs.dtype, np.int64)
+
+        p = np.poly1d(z.astype(np.float32))
+        assert_equal(p.coeffs.dtype, np.float32)
+
+        p = np.poly1d(z.astype(np.complex64))
+        assert_equal(p.coeffs.dtype, np.complex64)
+
     def test_poly_eq(self):
         p = np.poly1d([1, 2, 3])
         p2 = np.poly1d([1, 2, 4])


Code review comments:
diff_hunk: @@ -1236,7 +1236,7 @@ def __init__(self, c_or_r, r=False, variable=None):
             raise ValueError("Polynomial must be 1d only.")
         c_or_r = trim_zeros(c_or_r, trim='f')
         if len(c_or_r) == 0:
-            c_or_r = NX.array([0.])
+            c_or_r = NX.array([0.], dtype=c_or_r.dtype)
path: numpy/lib/polynomial.py
text: @user1:
May as well make this an integer literal now too:
```suggestion
            c_or_r = NX.array([0], dtype=c_or_r.dtype)
```