poly1d binary operations return inconsistent type with array input

poly1d from the old-style polynomial API returns objects of different types from binary operations between poly1d objects and arrays depending on the ordering of the operations:

>>> p = np.poly1d([1, 2, 3])
>>> a = np.array([1, 2, 3])
>>> p + a
poly1d([2, 4, 6])
>>> a + p
array([2, 4, 6])

Note that this issue only pertains to the old poly1d class. The recommended polynomial API from the np.polynomial package does not have this problem.

A fix for this was attempted in #15756, but subsequent discussion in that thread made clear that it was an inadequate solution as it introduced silent, backwards-incompatible behavior changes.

Numpy/Python version information:

1.19.0.dev0+965b41d 3.8.2 (default, Feb 26 2020, 22:21:03)
[GCC 9.2.1 20200130]

Start with the old-style poly1d class and reproduce the two binary-operation examples using the reported Python and NumPy versions. Review the discussion and attempted fix in #15756 before deciding on the intended behavior, with done defined as consistent operations that avoid silent backwards-incompatible changes.