BUG: Make Polynomial evaluation adhere to nep 50 (numpy#26550)

* Make sure evaluation of a Polynomial respects nep50 for scalar input:
```
import numpy as np

p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32))
print(type(p(2))) # np.float64
# correct, since the domain argument is the default which maps to [-1., 1.]

w=np.array([-1,1], dtype=np.float32)
p=np.polynomial.Polynomial(np.array([1, 2], dtype=np.float32), domain=w, window=w)
print(type(p(2))) # np.float32 (was float64 on main)
```

* Update documentation of the various polynomial classes for the updated domain and window (was changed in numpy#24568)
* Not addressed:
```
import numpy as np
arr = np.polydiv(1, np.float32(1))
arr.dtype # float64
```
The input here are polynomial coefficients, which are really an array and not a scalar. So the output type seems correct, even though `1` looks like a scalar input.

Author: eendebakpt

benchmarks/benchmarks/bench_polynomial.py

@@ -0,0 +1,29 @@
+from .common import Benchmark
+
+import numpy as np
+
+
+class Polynomial(Benchmark):
+
+    def setup(self):
+        self.polynomial_degree2 = np.polynomial.Polynomial(np.array([1, 2]))
+        self.array3 = np.linspace(0, 1, 3)
+        self.array1000 = np.linspace(0, 1, 10_000)
+        self.float64 = np.float64(1.0)
+
+    def time_polynomial_evaluation_scalar(self):
+        self.polynomial_degree2(self.float64)
+
+    def time_polynomial_evaluation_python_float(self):
+        self.polynomial_degree2(1.0)
+
+    def time_polynomial_evaluation_array_3(self):
+        self.polynomial_degree2(self.array3)
+
+    def time_polynomial_evaluation_array_1000(self):
+        self.polynomial_degree2(self.array1000)
+        
+    def time_polynomial_addition(self):
+        _ = self.polynomial_degree2 + self.polynomial_degree2
+        
+

numpy/polynomial/chebyshev.py

@@ -2011,9 +2011,9 @@ class Chebyshev(ABCPolyBase):
     domain : (2,) array_like, optional
         Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
         to the interval ``[window[0], window[1]]`` by shifting and scaling.
-        The default value is [-1, 1].
+        The default value is [-1., 1.].
     window : (2,) array_like, optional
-        Window, see `domain` for its use. The default value is [-1, 1].
+        Window, see `domain` for its use. The default value is [-1., 1.].
 
         .. versionadded:: 1.6.0
     symbol : str, optional

numpy/polynomial/hermite.py

@@ -1756,9 +1756,9 @@ class Hermite(ABCPolyBase):
     domain : (2,) array_like, optional
         Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
         to the interval ``[window[0], window[1]]`` by shifting and scaling.
-        The default value is [-1, 1].
+        The default value is [-1., 1.].
     window : (2,) array_like, optional
-        Window, see `domain` for its use. The default value is [-1, 1].
+        Window, see `domain` for its use. The default value is [-1., 1.].
 
         .. versionadded:: 1.6.0
     symbol : str, optional

numpy/polynomial/hermite_e.py

@@ -1667,9 +1667,9 @@ class HermiteE(ABCPolyBase):
     domain : (2,) array_like, optional
         Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
         to the interval ``[window[0], window[1]]`` by shifting and scaling.
-        The default value is [-1, 1].
+        The default value is [-1., 1.].
     window : (2,) array_like, optional
-        Window, see `domain` for its use. The default value is [-1, 1].
+        Window, see `domain` for its use. The default value is [-1., 1.].
 
         .. versionadded:: 1.6.0
     symbol : str, optional

numpy/polynomial/laguerre.py

@@ -1688,9 +1688,9 @@ class Laguerre(ABCPolyBase):
     domain : (2,) array_like, optional
         Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
         to the interval ``[window[0], window[1]]`` by shifting and scaling.
-        The default value is [0, 1].
+        The default value is [0., 1.].
     window : (2,) array_like, optional
-        Window, see `domain` for its use. The default value is [0, 1].
+        Window, see `domain` for its use. The default value is [0., 1.].
 
         .. versionadded:: 1.6.0
     symbol : str, optional

numpy/polynomial/legendre.py

@@ -1632,9 +1632,9 @@ class Legendre(ABCPolyBase):
     domain : (2,) array_like, optional
         Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
         to the interval ``[window[0], window[1]]`` by shifting and scaling.
-        The default value is [-1, 1].
+        The default value is [-1., 1.].
     window : (2,) array_like, optional
-        Window, see `domain` for its use. The default value is [-1, 1].
+        Window, see `domain` for its use. The default value is [-1., 1.].
 
         .. versionadded:: 1.6.0
     symbol : str, optional

numpy/polynomial/polynomial.py

@@ -1601,9 +1601,9 @@ class Polynomial(ABCPolyBase):
     domain : (2,) array_like, optional
         Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
         to the interval ``[window[0], window[1]]`` by shifting and scaling.
-        The default value is [-1, 1].
+        The default value is [-1., 1.].
     window : (2,) array_like, optional
-        Window, see `domain` for its use. The default value is [-1, 1].
+        Window, see `domain` for its use. The default value is [-1., 1.].
 
         .. versionadded:: 1.6.0
     symbol : str, optional

numpy/polynomial/polyutils.py

@@ -346,7 +346,8 @@ def mapdomain(x, old, new):
     array([-1.0+1.j , -0.6+0.6j, -0.2+0.2j,  0.2-0.2j,  0.6-0.6j,  1.0-1.j ]) # may vary
 
     """
-    x = np.asanyarray(x)
+    if type(x) not in (int, float, complex) and not isinstance(x, np.generic):
+        x = np.asanyarray(x)
     off, scl = mapparms(old, new)
     return off + scl*x
 

numpy/polynomial/tests/test_polynomial.py

@@ -627,3 +627,12 @@ def test_polyline(self):
 
     def test_polyline_zero(self):
         assert_equal(poly.polyline(3, 0), [3])
+
+    def test_result_type(self):
+        w = np.array([-1, 1], dtype=np.float32)
+        p = np.polynomial.Polynomial(w, domain=w, window=w)
+        v = p(2)
+        assert_equal(v.dtype, np.float32)
+
+        arr = np.polydiv(1, np.float32(1))
+        assert_equal(arr[0].dtype, np.float64)