MAINT: Unify polynomial division functions

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

Note that while chebdiv and polydiv could be implemented in terms of this function, their current implementations are more optimal and exploit the properties of a multiplication by a basis polynomial.

7 files changed, 42 insertions, 80 deletions

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index fd05280e9..0eef90177 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -793,6 +793,7 @@ def chebdiv(c1, c2):
     if c2[-1] == 0:
         raise ZeroDivisionError()
 
+    # note: this is more efficient than `pu._div(chebmul, c1, c2)`
     lc1 = len(c1)
     lc2 = len(c2)
     if lc1 < lc2:
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 605cb29ad..3767a80fc 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -550,26 +550,7 @@ def hermdiv(c1, c2):
     (array([1., 2., 3.]), array([1., 1.]))
 
     """
-    # c1, c2 are trimmed copies
-    [c1, c2] = pu.as_series([c1, c2])
-    if c2[-1] == 0:
-        raise ZeroDivisionError()
-
-    lc1 = len(c1)
-    lc2 = len(c2)
-    if lc1 < lc2:
-        return c1[:1]*0, c1
-    elif lc2 == 1:
-        return c1/c2[-1], c1[:1]*0
-    else:
-        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
-        rem = c1
-        for i in range(lc1 - lc2, - 1, -1):
-            p = hermmul([0]*i + [1], c2)
-            q = rem[-1]/p[-1]
-            rem = rem[:-1] - q*p[:-1]
-            quo[i] = q
-        return quo, pu.trimseq(rem)
+    return pu._div(hermmul, c1, c2)
 
 
 def hermpow(c, pow, maxpower=16):
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index b28881013..228396457 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -545,26 +545,7 @@ def hermediv(c1, c2):
     (array([1., 2., 3.]), array([1., 2.]))
 
     """
-    # c1, c2 are trimmed copies
-    [c1, c2] = pu.as_series([c1, c2])
-    if c2[-1] == 0:
-        raise ZeroDivisionError()
-
-    lc1 = len(c1)
-    lc2 = len(c2)
-    if lc1 < lc2:
-        return c1[:1]*0, c1
-    elif lc2 == 1:
-        return c1/c2[-1], c1[:1]*0
-    else:
-        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
-        rem = c1
-        for i in range(lc1 - lc2, - 1, -1):
-            p = hermemul([0]*i + [1], c2)
-            q = rem[-1]/p[-1]
-            rem = rem[:-1] - q*p[:-1]
-            quo[i] = q
-        return quo, pu.trimseq(rem)
+    return pu._div(hermemul, c1, c2)
 
 
 def hermepow(c, pow, maxpower=16):
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 575c5b2bc..dec469c17 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -547,26 +547,7 @@ def lagdiv(c1, c2):
     (array([1., 2., 3.]), array([1., 1.]))
 
     """
-    # c1, c2 are trimmed copies
-    [c1, c2] = pu.as_series([c1, c2])
-    if c2[-1] == 0:
-        raise ZeroDivisionError()
-
-    lc1 = len(c1)
-    lc2 = len(c2)
-    if lc1 < lc2:
-        return c1[:1]*0, c1
-    elif lc2 == 1:
-        return c1/c2[-1], c1[:1]*0
-    else:
-        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
-        rem = c1
-        for i in range(lc1 - lc2, - 1, -1):
-            p = lagmul([0]*i + [1], c2)
-            q = rem[-1]/p[-1]
-            rem = rem[:-1] - q*p[:-1]
-            quo[i] = q
-        return quo, pu.trimseq(rem)
+    return pu._div(lagmul, c1, c2)
 
 
 def lagpow(c, pow, maxpower=16):
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 6cd4360da..5f8ce53a3 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -590,26 +590,7 @@ def legdiv(c1, c2):
     (array([-0.07407407,  1.66666667]), array([-1.03703704, -2.51851852])) # may vary
 
     """
-    # c1, c2 are trimmed copies
-    [c1, c2] = pu.as_series([c1, c2])
-    if c2[-1] == 0:
-        raise ZeroDivisionError()
-
-    lc1 = len(c1)
-    lc2 = len(c2)
-    if lc1 < lc2:
-        return c1[:1]*0, c1
-    elif lc2 == 1:
-        return c1/c2[-1], c1[:1]*0
-    else:
-        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
-        rem = c1
-        for i in range(lc1 - lc2, - 1, -1):
-            p = legmul([0]*i + [1], c2)
-            q = rem[-1]/p[-1]
-            rem = rem[:-1] - q*p[:-1]
-            quo[i] = q
-        return quo, pu.trimseq(rem)
+    return pu._div(legmul, c1, c2)
 
 
 def legpow(c, pow, maxpower=16):
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 8c6b604d5..f63d9dd74 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -397,6 +397,7 @@ def polydiv(c1, c2):
     if c2[-1] == 0:
         raise ZeroDivisionError()
 
+    # note: this is more efficient than `pu._div(polymul, c1, c2)`
     lc1 = len(c1)
     lc2 = len(c2)
     if lc1 < lc2:
diff --git a/numpy/polynomial/polyutils.py b/numpy/polynomial/polyutils.py
index db1cb2841..1d5390984 100644
--- a/numpy/polynomial/polyutils.py
+++ b/numpy/polynomial/polyutils.py
@@ -537,3 +537,39 @@ def _gridnd(val_f, c, *args):
     for xi in args:
         c = val_f(xi, c)
     return c
+
+
+def _div(mul_f, c1, c2):
+    """
+    Helper function used to implement the ``<type>div`` functions.
+
+    Implementation uses repeated subtraction of c2 multiplied by the nth basis.
+    For some polynomial types, a more efficient approach may be possible.
+
+    Parameters
+    ----------
+    mul_f : function(array_like, array_like) -> array_like
+        The ``<type>mul`` function, such as ``polymul``
+    c1, c2 :
+        See the ``<type>div`` functions for more detail
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = as_series([c1, c2])
+    if c2[-1] == 0:
+        raise ZeroDivisionError()
+
+    lc1 = len(c1)
+    lc2 = len(c2)
+    if lc1 < lc2:
+        return c1[:1]*0, c1
+    elif lc2 == 1:
+        return c1/c2[-1], c1[:1]*0
+    else:
+        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
+        rem = c1
+        for i in range(lc1 - lc2, - 1, -1):
+            p = mul_f([0]*i + [1], c2)
+            q = rem[-1]/p[-1]
+            rem = rem[:-1] - q*p[:-1]
+            quo[i] = q
+        return quo, trimseq(rem)