MAINT: Unify polynomial power functions

These power functions are all the same - the algorithm used does not care about the basis.

`polypow` and `chebpow` have some optimizations in their versions, which this maintains

7 files changed, 43 insertions, 90 deletions

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 475582fdc..836f47363 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -826,6 +826,9 @@ def chebpow(c, pow, maxpower=16):
     array([15.5, 22. , 16. , ..., 12.5, 12. ,  8. ])
 
     """
+    # note: this is more efficient than `pu._pow(chebmul, c1, c2)`, as it
+    # avoids converting between z and c series repeatedly
+
     # c is a trimmed copy
     [c] = pu.as_series([c])
     power = int(pow)
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index e925c4304..1f9f07d29 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -570,24 +570,7 @@ def hermpow(c, pow, maxpower=16):
     array([81.,  52.,  82.,  12.,   9.])
 
     """
-    # c is a trimmed copy
-    [c] = pu.as_series([c])
-    power = int(pow)
-    if power != pow or power < 0:
-        raise ValueError("Power must be a non-negative integer.")
-    elif maxpower is not None and power > maxpower:
-        raise ValueError("Power is too large")
-    elif power == 0:
-        return np.array([1], dtype=c.dtype)
-    elif power == 1:
-        return c
-    else:
-        # This can be made more efficient by using powers of two
-        # in the usual way.
-        prd = c
-        for i in range(2, power + 1):
-            prd = hermmul(prd, c)
-        return prd
+    return pu._pow(hermmul, c, pow, maxpower)
 
 
 def hermder(c, m=1, scl=1, axis=0):
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index 8ae8c7f95..48e5eab20 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -565,24 +565,7 @@ def hermepow(c, pow, maxpower=16):
     array([23.,  28.,  46.,  12.,   9.])
 
     """
-    # c is a trimmed copy
-    [c] = pu.as_series([c])
-    power = int(pow)
-    if power != pow or power < 0:
-        raise ValueError("Power must be a non-negative integer.")
-    elif maxpower is not None and power > maxpower:
-        raise ValueError("Power is too large")
-    elif power == 0:
-        return np.array([1], dtype=c.dtype)
-    elif power == 1:
-        return c
-    else:
-        # This can be made more efficient by using powers of two
-        # in the usual way.
-        prd = c
-        for i in range(2, power + 1):
-            prd = hermemul(prd, c)
-        return prd
+    return pu._pow(hermemul, c, pow, maxpower)
 
 
 def hermeder(c, m=1, scl=1, axis=0):
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 89e053ffa..41d15142a 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -567,24 +567,7 @@ def lagpow(c, pow, maxpower=16):
     array([ 14., -16.,  56., -72.,  54.])
 
     """
-    # c is a trimmed copy
-    [c] = pu.as_series([c])
-    power = int(pow)
-    if power != pow or power < 0:
-        raise ValueError("Power must be a non-negative integer.")
-    elif maxpower is not None and power > maxpower:
-        raise ValueError("Power is too large")
-    elif power == 0:
-        return np.array([1], dtype=c.dtype)
-    elif power == 1:
-        return c
-    else:
-        # This can be made more efficient by using powers of two
-        # in the usual way.
-        prd = c
-        for i in range(2, power + 1):
-            prd = lagmul(prd, c)
-        return prd
+    return pu._pow(lagmul, c, pow, maxpower)
 
 
 def lagder(c, m=1, scl=1, axis=0):
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 4b56f82f4..d56e2ae2c 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -607,24 +607,7 @@ def legpow(c, pow, maxpower=16):
     --------
 
     """
-    # c is a trimmed copy
-    [c] = pu.as_series([c])
-    power = int(pow)
-    if power != pow or power < 0:
-        raise ValueError("Power must be a non-negative integer.")
-    elif maxpower is not None and power > maxpower:
-        raise ValueError("Power is too large")
-    elif power == 0:
-        return np.array([1], dtype=c.dtype)
-    elif power == 1:
-        return c
-    else:
-        # This can be made more efficient by using powers of two
-        # in the usual way.
-        prd = c
-        for i in range(2, power + 1):
-            prd = legmul(prd, c)
-        return prd
+    return pu._pow(legmul, c, pow, maxpower)
 
 
 def legder(c, m=1, scl=1, axis=0):
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 44e1cbb3f..e3eb7ec5a 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -434,24 +434,9 @@ def polypow(c, pow, maxpower=None):
     array([ 1., 4., 10., 12., 9.])
 
     """
-    # c is a trimmed copy
-    [c] = pu.as_series([c])
-    power = int(pow)
-    if power != pow or power < 0:
-        raise ValueError("Power must be a non-negative integer.")
-    elif maxpower is not None and power > maxpower:
-        raise ValueError("Power is too large")
-    elif power == 0:
-        return np.array([1], dtype=c.dtype)
-    elif power == 1:
-        return c
-    else:
-        # This can be made more efficient by using powers of two
-        # in the usual way.
-        prd = c
-        for i in range(2, power + 1):
-            prd = np.convolve(prd, c)
-        return prd
+    # note: this is more efficient than `pu._pow(polymul, c1, c2)`, as it
+    # avoids calling `as_series` repeatedly
+    return pu._pow(np.convolve, c, pow, maxpower)
 
 
 def polyder(c, m=1, scl=1, axis=0):
diff --git a/numpy/polynomial/polyutils.py b/numpy/polynomial/polyutils.py
index 1a73f47d8..5128e35e0 100644
--- a/numpy/polynomial/polyutils.py
+++ b/numpy/polynomial/polyutils.py
@@ -690,6 +690,39 @@ def _fit(vander_f, x, y, deg, rcond=None, full=False, w=None):
         return c
 
 
+def _pow(mul_f, c, pow, maxpower):
+    """
+    Helper function used to implement the ``<type>pow`` functions.
+
+    Parameters
+    ----------
+    vander_f : function(array_like, int) -> ndarray
+        The 1d vander function, such as ``polyvander``
+    pow, maxpower :
+        See the ``<type>pow`` functions for more detail
+    mul_f : function(array_like, array_like) -> ndarray
+        The ``<type>mul`` function, such as ``polymul``
+    """
+    # c is a trimmed copy
+    [c] = as_series([c])
+    power = int(pow)
+    if power != pow or power < 0:
+        raise ValueError("Power must be a non-negative integer.")
+    elif maxpower is not None and power > maxpower:
+        raise ValueError("Power is too large")
+    elif power == 0:
+        return np.array([1], dtype=c.dtype)
+    elif power == 1:
+        return c
+    else:
+        # This can be made more efficient by using powers of two
+        # in the usual way.
+        prd = c
+        for i in range(2, power + 1):
+            prd = mul_f(prd, c)
+        return prd
+
+
 def _deprecate_as_int(x, desc):
     """
     Like `operator.index`, but emits a deprecation warning when passed a float