Fix polynomial division

2 files changed, 27 insertions, 20 deletions

diff --git a/scipy/base/polynomial.py b/scipy/base/polynomial.py
index 272f20f80..eb0856676 100644
--- a/scipy/base/polynomial.py
+++ b/scipy/base/polynomial.py
@@ -240,7 +240,8 @@ def polyadd(a1, a2):
     """Adds two polynomials represented as sequences
     """
     truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
-    a1, a2 = map(atleast_1d, (a1, a2))
+    a1 = atleast_1d(a1)
+    a2 = atleast_1d(a2)
     diff = len(a2) - len(a1)
     if diff == 0:
         return a1 + a2
@@ -258,7 +259,8 @@ def polysub(a1, a2):
     """Subtracts two polynomials represented as sequences
     """
     truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
-    a1, a2 = map(atleast_1d, (a1, a2))
+    a1 = atleast_1d(a1)
+    a2 = atleast_1d(a2)
     diff = len(a2) - len(a1)
     if diff == 0:
         return a1 - a2
@@ -301,26 +303,29 @@ def deconvolve(signal, divisor):
         rem = num - NX.convolve(den, quot, mode='full')
     return quot, rem
 
-def polydiv(a1, a2):
-    """Computes q and r polynomials so that a1(s) = q(s)*a2(s) + r(s)
+def polydiv(u, v):
+    """Computes q and r polynomials so that u(s) = q(s)*v(s) + r(s)
+    and deg r < deg v.
     """
-    truepoly = (isinstance(a1, poly1d) or isinstance(a2, poly1d))
-    monDivisor = NX.asarray(a2) / a2[0]
-    dividend = NX.asarray(a1) / a2[0]
-    q = []
-    r = dividend
-    while len(r) >= len(monDivisor):
-        q.append(r[0])
-        r = polysub(r, polymul(q, monDivisor))[1:]
-    q = NX.asarray(q)
+    truepoly = (isinstance(u, poly1d) or isinstance(u, poly1d))
+    u = atleast_1d(u)
+    v = atleast_1d(v)
+    m = len(u) - 1
+    n = len(v) - 1
+    scale = 1. / v[0]
+    q = NX.zeros((m-n+1,), float)
+    r = u.copy()
+    for k in range(0, m-n+1):
+        d = scale * r[k]
+        q[k] = d
+        r[k:k+n+1] -= d*v
     while NX.allclose(r[0], 0, rtol=1e-14) and (r.shape[-1] > 1):
         r = r[1:]
-    r *= a2[0]
     if truepoly:
-        q, r = map(poly1d, (q, r))
+        q = poly1d(q)
+        r = poly1d(r)
     return q, r
 
-
 _poly_mat = re.compile(r"[*][*]([0-9]*)")
 def _raise_power(astr, wrap=70):
     n = 0
@@ -498,14 +503,14 @@ class poly1d(object):
             return poly1d(self.coeffs/other)
         else:
             other = poly1d(other)
-            return map(poly1d, polydiv(self.coeffs, other.coeffs))
+            return polydiv(self, other)
 
     def __rdiv__(self, other):
         if isscalar(other):
             return poly1d(other/self.coeffs)
         else:
             other = poly1d(other)
-            return map(poly1d, polydiv(other.coeffs, self.coeffs))
+            return polydiv(other, self)
 
     def __setattr__(self, key, val):
         raise ValueError, "Attributes cannot be changed this way."
diff --git a/scipy/base/tests/test_polynomial.py b/scipy/base/tests/test_polynomial.py
index 928306e32..421368e67 100644
--- a/scipy/base/tests/test_polynomial.py
+++ b/scipy/base/tests/test_polynomial.py
@@ -1,6 +1,6 @@
 """
 >>> import scipy.base as nx
->>> from scipy.base.polynomial import poly1d
+>>> from scipy.base.polynomial import poly1d, polydiv
 
 >>> p = poly1d([1.,2,3])
 >>> p
@@ -27,7 +27,7 @@ poly1d([ 3.,  2.,  1.])
 >>> p * q
 poly1d([  3.,   8.,  14.,   8.,   3.])
 >>> p / q
-[poly1d([ 0.33333333]), poly1d([ 1.33333333,  2.66666667])]
+(poly1d([ 0.33333333]), poly1d([ 1.33333333,  2.66666667]))
 >>> p + q
 poly1d([ 4.,  4.,  4.])
 >>> p - q
@@ -69,6 +69,8 @@ poly1d([ 2.])
         2
 1 lambda + 2 lambda + 3
 
+>>> polydiv(poly1d([1,0,-1]), poly1d([1,1]))
+(poly1d([ 1., -1.]), poly1d([ 0.]))
 """
 
 from scipy.test.testing import *