ENH: Add companion matrix functions.

The new companion matrices are related to the old by a
similarity transformation that makes them better conditioned
for root finding. In particular, the companion matrices for
the orthogonal polynomials are symmetric when the zeros of a
single polynomial term is wanted. This produces better zeros
for use in Gauss quadrature.

1 files changed, 50 insertions, 13 deletions

diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 345d76233..1d8c1dbd3 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -1246,6 +1246,36 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
         return c
 
 
+def polycompanion(cs):
+    """Return the companion matrix of cs.
+
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of series coefficients ordered from low to high degree.
+
+    Returns
+    -------
+    mat : ndarray
+        Scaled companion matrix of dimensions (deg, deg).
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) < 2 :
+        raise ValueError('Series must have maximum degree of at least 1.')
+    if len(cs) == 2:
+        return np.array(-cs[0]/cs[1])
+
+    n = len(cs) - 1
+    mat = np.zeros((n, n), dtype=cs.dtype)
+    bot = mat.reshape(-1)[n::n+1]
+    bot[...] = 1
+    mat[:,-1] -= cs[:-1]/cs[-1]
+    return mat
+
+
 def polyroots(cs):
     """
     Compute the roots of a polynomial.
@@ -1270,32 +1300,39 @@ def polyroots(cs):
     --------
     chebroots
 
+    Notes
+    -----
+    The root estimates are obtained as the eigenvalues of the companion
+    matrix, Roots far from the origin of the complex plane may have large
+    errors due to the numerical instability of the power series for such
+    values. Roots with multiplicity greater than 1 will also show larger
+    errors as the value of the series near such points is relatively
+    insensitive to errors in the roots. Isolated roots near the origin can
+    be improved by a few iterations of Newton's method.
+
     Examples
     --------
-    >>> import numpy.polynomial as P
-    >>> P.polyroots(P.polyfromroots((-1,0,1)))
+    >>> import numpy.polynomial.polynomial as poly
+    >>> poly.polyroots(poly.polyfromroots((-1,0,1)))
     array([-1.,  0.,  1.])
-    >>> P.polyroots(P.polyfromroots((-1,0,1))).dtype
+    >>> poly.polyroots(poly.polyfromroots((-1,0,1))).dtype
     dtype('float64')
     >>> j = complex(0,1)
-    >>> P.polyroots(P.polyfromroots((-j,0,j)))
+    >>> poly.polyroots(poly.polyfromroots((-j,0,j)))
     array([  0.00000000e+00+0.j,   0.00000000e+00+1.j,   2.77555756e-17-1.j])
 
     """
     # cs is a trimmed copy
     [cs] = pu.as_series([cs])
-    if len(cs) <= 1 :
+    if len(cs) < 2:
         return np.array([], dtype=cs.dtype)
-    if len(cs) == 2 :
+    if len(cs) == 2:
         return np.array([-cs[0]/cs[1]])
 
-    n = len(cs) - 1
-    cmat = np.zeros((n,n), dtype=cs.dtype)
-    cmat.flat[n::n+1] = 1
-    cmat[:,-1] -= cs[:-1]/cs[-1]
-    roots = la.eigvals(cmat)
-    roots.sort()
-    return roots
+    m = polycompanion(cs)
+    r = la.eigvals(m)
+    r.sort()
+    return r
 
 
 #