ENH: Improve the computation of polynomials from roots.

The original method was overly sensitive to roundoff. Of the two
approaches considered, gauss integration or binary subdivision of
the roots, the latter is more compatible with using other number
representations such as mpmath. No method is going to be suitable
for large numbers of arbitrary zeros but the current method is a
significant improvement.

6 files changed, 66 insertions, 24 deletions

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 022dccfa7..ab4c96a54 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -540,10 +540,17 @@ def chebfromroots(roots) :
         return np.ones(1)
     else :
         [roots] = pu.as_series([roots], trim=False)
-        prd = np.array([1], dtype=roots.dtype)
-        for r in roots:
-            prd = chebsub(chebmulx(prd), r*prd)
-        return prd
+        roots.sort()
+        n = len(roots)
+        p = [chebline(-r, 1) for r in roots]
+        while n > 1:
+            m = n//2
+            tmp = [chebmul(p[i], p[i+m]) for i in range(m)]
+            if n%2:
+                tmp[0] = chebmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def chebadd(c1, c2):
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 627232d7c..32b50a866 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -287,10 +287,17 @@ def hermfromroots(roots) :
         return np.ones(1)
     else :
         [roots] = pu.as_series([roots], trim=False)
-        prd = np.array([1], dtype=roots.dtype)
-        for r in roots:
-            prd = hermsub(hermmulx(prd), r*prd)
-        return prd
+        roots.sort()
+        n = len(roots)
+        p = [hermline(-r, 1) for r in roots]
+        while n > 1:
+            m = n//2
+            tmp = [hermmul(p[i], p[i+m]) for i in range(m)]
+            if n%2:
+                tmp[0] = hermmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def hermadd(c1, c2):
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index daed05008..0511307a4 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -287,10 +287,17 @@ def hermefromroots(roots) :
         return np.ones(1)
     else :
         [roots] = pu.as_series([roots], trim=False)
-        prd = np.array([1], dtype=roots.dtype)
-        for r in roots:
-            prd = hermesub(hermemulx(prd), r*prd)
-        return prd
+        roots.sort()
+        n = len(roots)
+        p = [hermeline(-r, 1) for r in roots]
+        while n > 1:
+            m = n//2
+            tmp = [hermemul(p[i], p[i+m]) for i in range(m)]
+            if n%2:
+                tmp[0] = hermemul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def hermeadd(c1, c2):
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index fc9afec00..866aab5c4 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -283,10 +283,17 @@ def lagfromroots(roots) :
         return np.ones(1)
     else :
         [roots] = pu.as_series([roots], trim=False)
-        prd = np.array([1], dtype=roots.dtype)
-        for r in roots:
-            prd = lagsub(lagmulx(prd), r*prd)
-        return prd
+        roots.sort()
+        n = len(roots)
+        p = [lagline(-r, 1) for r in roots]
+        while n > 1:
+            m = n//2
+            tmp = [lagmul(p[i], p[i+m]) for i in range(m)]
+            if n%2:
+                tmp[0] = lagmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def lagadd(c1, c2):
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 8fd64985b..b76e678c2 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -314,10 +314,17 @@ def legfromroots(roots) :
         return np.ones(1)
     else :
         [roots] = pu.as_series([roots], trim=False)
-        prd = np.array([1], dtype=roots.dtype)
-        for r in roots:
-            prd = legsub(legmulx(prd), r*prd)
-        return prd
+        roots.sort()
+        n = len(roots)
+        p = [legline(-r, 1) for r in roots]
+        while n > 1:
+            m = n//2
+            tmp = [legmul(p[i], p[i+m]) for i in range(m)]
+            if n%2:
+                tmp[0] = legmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def legadd(c1, c2):
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 207511216..c8894c68e 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -186,10 +186,17 @@ def polyfromroots(roots) :
         return np.ones(1)
     else :
         [roots] = pu.as_series([roots], trim=False)
-        prd = np.array([1], dtype=roots.dtype)
-        for r in roots:
-            prd = polysub(polymulx(prd), r*prd)
-        return prd
+        roots.sort()
+        n = len(roots)
+        p = [polyline(-r, 1) for r in roots]
+        while n > 1:
+            m = n//2
+            tmp = [polymul(p[i], p[i+m]) for i in range(m)]
+            if n%2:
+                tmp[0] = polymul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def polyadd(c1, c2):