Merge pull request #263 from charris/backport-poly

Backport numpy.poly bug fixes and documentation improvements.

18 files changed, 1594 insertions, 317 deletions

diff --git a/doc/source/reference/routines.polynomials.chebyshev.rst b/doc/source/reference/routines.polynomials.chebyshev.rst
new file mode 100644
index 000000000..60c816f03
--- /dev/null
+++ b/doc/source/reference/routines.polynomials.chebyshev.rst
@@ -0,0 +1,92 @@
+Chebyshev Module (:mod:`numpy.polynomial.chebyshev`)
+====================================================
+
+.. versionadded:: 1.4.0
+
+.. currentmodule:: numpy.polynomial.chebyshev
+
+This module provides a number of objects (mostly functions) useful for
+dealing with Chebyshev series, including a `Chebyshev` class that
+encapsulates the usual arithmetic operations.  (General information
+on how this module represents and works with such polynomials is in the
+docstring for its "parent" sub-package, `numpy.polynomial`).
+
+Chebyshev Class
+---------------
+
+.. autosummary::
+   :toctree: generated/
+
+   Chebyshev
+
+Basics
+------
+
+.. autosummary::
+   :toctree: generated/
+
+   chebval
+   chebval2d
+   chebval3d
+   chebgrid2d
+   chebgrid3d
+   chebroots
+   chebfromroots
+
+Fitting
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   chebfit
+   chebvander
+   chebvander2d
+   chebvander3d
+
+Calculus
+--------
+
+.. autosummary::
+   :toctree: generated/
+
+   chebder
+   chebint
+
+Algebra
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   chebadd
+   chebsub
+   chebmul
+   chebmulx
+   chebdiv
+   chebpow
+
+Quadrature
+----------
+
+.. autosummary::
+   :toctree: generated/
+
+   chebgauss
+   chebweight
+
+Miscellaneous
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   chebcompanion
+   chebdomain
+   chebzero
+   chebone
+   chebx
+   chebtrim
+   chebline
+   cheb2poly
+   poly2cheb
diff --git a/doc/source/reference/routines.polynomials.classes.rst b/doc/source/reference/routines.polynomials.classes.rst
new file mode 100644
index 000000000..9139b802e
--- /dev/null
+++ b/doc/source/reference/routines.polynomials.classes.rst
@@ -0,0 +1,316 @@
+Using the Convenience Classes
+=============================
+
+The convenience classes provided by the polynomial package are:
+
+         ============    ================
+         Name            Provides
+         ============    ================
+         Polynomial      Power series
+         Chebyshev       Chebyshev series
+         Legendre        Legendre series
+         Laguerre        Laguerre series
+         Hermite         Hermite series
+         HermiteE        HermiteE series
+         ============    ================
+
+The series in this context are finite sums of the corresponding polynomial
+basis functions multiplied by coefficients. For instance, a power series
+looks like
+
+.. math:: p(x) = 1 + 2x + 3x^2
+
+and has coefficients :math:`[1, 2, 3]`. The Chebyshev series with the
+same coefficients looks like
+
+
+.. math:: p(x) = 1 T_0(x) + 2 T_1(x) + 3 T_2(x)
+
+and more generally
+
+.. math:: p(x) = \sum_{i=0}^n c_i T_i(x)
+
+where in this case the :math:`T_n` are the Chebyshev functions of degree
+`n`, but could just as easily be the basis functions of any of the other
+classes. The convention for all the classes is that the coefficient
+c[i] goes with the basis function of degree i.
+
+All of the classes have the same methods, and especially they implement the
+Python numeric operators +, -, \*, //, %, divmod, \*\*, ==,
+and !=. The last two can be a bit problematic due to floating point
+roundoff errors. We now give a quick demonstration of the various
+operations using Numpy version 1.7.0.
+
+Basics
+------
+
+First we need a polynomial class and a polynomial instance to play with.
+The classes can be imported directly from the polynomial package or from
+the module of the relevant type. Here we import from the package and use
+the conventional Polynomial class because of its familiarity.::
+
+   >>> from numpy.polynomial import Polynomial as P
+   >>> p = P([1,2,3])
+   >>> p
+   Polynomial([ 1.,  2.,  3.], [-1.,  1.], [-1.,  1.])
+
+Note that there are three parts to the long version of the printout. The
+first is the coefficients, the second is the domain, and the third is the
+window::
+
+   >>> p.coef
+   array([ 1.,  2.,  3.])
+   >>> p.domain
+   array([-1.,  1.])
+   >>> p.window
+   array([-1.,  1.])
+
+Printing a polynomial yields a shorter form without the domain
+and window::
+
+   >>> print p
+   poly([ 1.  2.  3.])
+
+We will deal with the domain and window when we get to fitting, for the moment
+we ignore them and run through the basic algebraic and arithmetic operations.
+
+Addition and Subtraction::
+
+   >>> p + p
+   Polynomial([ 2.,  4.,  6.], [-1.,  1.], [-1.,  1.])
+   >>> p - p
+   Polynomial([ 0.], [-1.,  1.], [-1.,  1.])
+
+Multiplication::
+
+   >>> p * p
+   Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])
+
+Powers::
+
+   >>> p**2
+   Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])
+
+Division:
+
+Floor_division, '//', is the division operator for the polynomial classes,
+polynomials are treated like integers in this regard. For Python versions <
+3.x the '/' operator maps to '//', as it does for Python, for later
+versions the '/' will only work for division by scalars. At some point it
+will be deprecated::
+
+   >>> p // P([-1, 1])
+   Polynomial([ 5.,  3.], [-1.,  1.], [-1.,  1.])
+
+Remainder::
+
+   >>> p % P([-1, 1])
+   Polynomial([ 6.], [-1.,  1.], [-1.,  1.])
+
+Divmod::
+
+   >>> quo, rem = divmod(p, P([-1, 1]))
+   >>> quo
+   Polynomial([ 5.,  3.], [-1.,  1.], [-1.,  1.])
+   >>> rem
+   Polynomial([ 6.], [-1.,  1.], [-1.,  1.])
+
+Evaluation::
+
+   >>> x = np.arange(5)
+   >>> p(x)
+   array([  1.,   6.,  17.,  34.,  57.])
+   >>> x = np.arange(6).reshape(3,2)
+   >>> p(x)
+   array([[  1.,   6.],
+          [ 17.,  34.],
+          [ 57.,  86.]])
+
+Substitution:
+
+Substitute a polynomial for x and expand the result. Here we substitute
+p in itself leading to a new polynomial of degree 4 after expansion. If
+the polynomials are regarded as functions this is composition of
+functions::
+
+   >>> p(p)
+   Polynomial([  6.,  16.,  36.,  36.,  27.], [-1.,  1.], [-1.,  1.])
+
+Roots::
+
+   >>> p.roots()
+   array([-0.33333333-0.47140452j, -0.33333333+0.47140452j])
+
+
+
+It isn't always convenient to explicitly use Polynomial instances, so
+tuples, lists, arrays, and scalars are automatically cast in the arithmetic
+operations::
+
+   >>> p + [1, 2, 3]
+   Polynomial([ 2.,  4.,  6.], [-1.,  1.], [-1.,  1.])
+   >>> [1, 2, 3] * p
+   Polynomial([  1.,   4.,  10.,  12.,   9.], [-1.,  1.], [-1.,  1.])
+   >>> p / 2
+   Polynomial([ 0.5,  1. ,  1.5], [-1.,  1.], [-1.,  1.])
+
+Polynomials that differ in domain, window, or class can't be mixed in
+arithmetic::
+
+    >>> from numpy.polynomial import Chebyshev as T
+    >>> p + P([1], domain=[0,1])
+    Traceback (most recent call last):
+      File "<stdin>", line 1, in <module>
+      File "<string>", line 213, in __add__
+    TypeError: Domains differ
+    >>> p + P([1], window=[0,1])
+    Traceback (most recent call last):
+      File "<stdin>", line 1, in <module>
+      File "<string>", line 215, in __add__
+    TypeError: Windows differ
+    >>> p + T([1])
+    Traceback (most recent call last):
+      File "<stdin>", line 1, in <module>
+      File "<string>", line 211, in __add__
+    TypeError: Polynomial types differ
+
+
+But different types can be used for substitution. In fact, this is how
+conversion of Polynomial classes among themselves is done for type, domain,
+and window casting::
+
+    >>> p(T([0, 1]))
+    Chebyshev([ 2.5,  2. ,  1.5], [-1.,  1.], [-1.,  1.])
+
+Which gives the polynomial 'p' in Chebyshev form. This works because
+:math:`T_1(x) = x` and substituting :math:`x` for :math:`x` doesn't change
+the original polynomial. However, all the multiplications and divisions
+will be done using Chebyshev series, hence the type of the result.
+
+Calculus
+--------
+
+Polynomial instances can be integrated and differentiated.::
+
+    >>> from numpy.polynomial import Polynomial as P
+    >>> p = P([2, 6])
+    >>> p.integ()
+    Polynomial([ 0.,  2.,  3.], [-1.,  1.], [-1.,  1.])
+    >>> p.integ(2)
+    Polynomial([ 0.,  0.,  1.,  1.], [-1.,  1.], [-1.,  1.])
+
+The first example integrates 'p' once, the second example integrates it
+twice. By default, the lower bound of the integration and the integration
+constant are 0, but both can be specified.::
+
+    >>> p.integ(lbnd=-1)
+    Polynomial([-1.,  2.,  3.], [-1.,  1.], [-1.,  1.])
+    >>> p.integ(lbnd=-1, k=1)
+    Polynomial([ 0.,  2.,  3.], [-1.,  1.], [-1.,  1.])
+
+In the first case the lower bound of the integration is set to -1 and the
+integration constant is 0. In the second the constant of integration is set
+to 1 as well. Differentiation is simpler since the only option is the
+number times the polynomial is differentiated::
+
+    >>> p = P([1, 2, 3])
+    >>> p.deriv(1)
+    Polynomial([ 2.,  6.], [-1.,  1.], [-1.,  1.])
+    >>> p.deriv(2)
+    Polynomial([ 6.], [-1.,  1.], [-1.,  1.])
+
+
+Other Polynomial Constructors
+-----------------------------
+
+Constructing polynomials by specifying coefficients is just one way of
+obtaining a polynomial instance, they may also be created by specifying
+their roots, by conversion from other polynomial types, and by least
+squares fits. Fitting is discussed in its own section, the other methods
+are demonstrated below.::
+
+    >>> from numpy.polynomial import Polynomial as P
+    >>> from numpy.polynomial import Chebyshev as T
+    >>> p = P.fromroots([1, 2, 3])
+    >>> p
+    Polynomial([ -6.,  11.,  -6.,   1.], [-1.,  1.], [-1.,  1.])
+    >>> p.convert(kind=T)
+    Chebyshev([ -9.  ,  11.75,  -3.  ,   0.25], [-1.,  1.], [-1.,  1.])
+
+The convert method can also convert domain and window::
+
+    >>> p.convert(kind=T, domain=[0, 1])
+    Chebyshev([-2.4375 ,  2.96875, -0.5625 ,  0.03125], [ 0.,  1.], [-1.,  1.])
+    >>> p.convert(kind=P, domain=[0, 1])
+    Polynomial([-1.875,  2.875, -1.125,  0.125], [ 0.,  1.], [-1.,  1.])
+
+Conversions between types can be useful, but it is *not* recommended
+for routine use. The loss of numerical precision in passing from a
+Chebyshev series of degree 50 to a Polynomial series of the same degree
+can make the results of numerical evaluation essentially random.
+
+Fitting
+-------
+
+Fitting is the reason that the `domain` and `window` attributes are part of
+the convenience classes. To illustrate the problem, the values of the Chebyshev
+polynomials up to degree 5 are plotted below.
+
+.. plot::
+
+    >>> import matplotlib.pyplot as plt
+    >>> from numpy.polynomial import Chebyshev as T
+    >>> x = np.linspace(-1, 1, 100)
+    >>> for i in range(6): ax = plt.plot(x, T.basis(i)(x), lw=2, label="T_%d"%i)
+    ...
+    >>> plt.legend(loc="upper left")
+    <matplotlib.legend.Legend object at 0x3b3ee10>
+    >>> plt.show()
+
+In the range -1 <= x <= 1 they are nice, equiripple functions lying between +/- 1.
+The same plots over the range -2 <= x <= 2 look very different:
+
+.. plot::
+
+    >>> import matplotlib.pyplot as plt
+    >>> from numpy.polynomial import Chebyshev as T
+    >>> x = np.linspace(-2, 2, 100)
+    >>> for i in range(6): ax = plt.plot(x, T.basis(i)(x), lw=2, label="T_%d"%i)
+    ...
+    >>> plt.legend(loc="lower right")
+    <matplotlib.legend.Legend object at 0x3b3ee10>
+    >>> plt.show()
+
+As can be seen, the "good" parts have shrunk to insignificance. In using
+Chebyshev polynomials for fitting we want to use the region where x is
+between -1 and 1 and that is what the 'window' specifies. However, it is
+unlikely that the data to be fit has all its data points in that interval,
+so we use 'domain' to specify the interval where the data points lie. When
+the fit is done, the domain is first mapped to the window by a linear
+transformation and the usual least squares fit is done using the mapped
+data points. The window and domain of the fit are part of the returned series
+and are automatically used when computing values, derivatives, and such. If
+they aren't specified in the call the fitting routine will use the default
+window and the smallest domain that holds all the data points. This is
+illustrated below for a fit to a noisy sin curve.
+
+.. plot::
+
+    >>> import numpy as np
+    >>> import matplotlib.pyplot as plt
+    >>> from numpy.polynomial import Chebyshev as T
+    >>> np.random.seed(11)
+    >>> x = np.linspace(0, 2*np.pi, 20)
+    >>> y = np.sin(x) + np.random.normal(scale=.1, size=x.shape)
+    >>> p = T.fit(x, y, 5)
+    >>> plt.plot(x, y, 'o')
+    [<matplotlib.lines.Line2D object at 0x2136c10>]
+    >>> xx, yy = p.linspace()
+    >>> plt.plot(xx, yy, lw=2)
+    [<matplotlib.lines.Line2D object at 0x1cf2890>]
+    >>> p.domain
+    array([ 0.        ,  6.28318531])
+    >>> p.window
+    array([-1.,  1.])
+    >>> plt.show()
+
diff --git a/doc/source/reference/routines.polynomials.hermite.rst b/doc/source/reference/routines.polynomials.hermite.rst
new file mode 100644
index 000000000..8ee72e97c
--- /dev/null
+++ b/doc/source/reference/routines.polynomials.hermite.rst
@@ -0,0 +1,92 @@
+Hermite Module, "Physicists'" (:mod:`numpy.polynomial.hermite`)
+===============================================================
+
+.. versionadded:: 1.6.0
+
+.. currentmodule:: numpy.polynomial.hermite
+
+This module provides a number of objects (mostly functions) useful for
+dealing with Hermite series, including a `Hermite` class that
+encapsulates the usual arithmetic operations.  (General information
+on how this module represents and works with such polynomials is in the
+docstring for its "parent" sub-package, `numpy.polynomial`).
+
+Hermite Class
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   Hermite
+
+Basics
+------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermval
+   hermval2d
+   hermval3d
+   hermgrid2d
+   hermgrid3d
+   hermroots
+   hermfromroots
+
+Fitting
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermfit
+   hermvander
+   hermvander2d
+   hermvander3d
+
+Calculus
+--------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermder
+   hermint
+
+Algebra
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermadd
+   hermsub
+   hermmul
+   hermmulx
+   hermdiv
+   hermpow
+
+Quadrature
+----------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermgauss
+   hermweight
+
+Miscellaneous
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermcompanion
+   hermdomain
+   hermzero
+   hermone
+   hermx
+   hermtrim
+   hermline
+   herm2poly
+   poly2herm
diff --git a/doc/source/reference/routines.polynomials.hermite_e.rst b/doc/source/reference/routines.polynomials.hermite_e.rst
new file mode 100644
index 000000000..33a15bb44
--- /dev/null
+++ b/doc/source/reference/routines.polynomials.hermite_e.rst
@@ -0,0 +1,92 @@
+HermiteE Module, "Probabilists'" (:mod:`numpy.polynomial.hermite_e`)
+====================================================================
+
+.. versionadded:: 1.6.0
+
+.. currentmodule:: numpy.polynomial.hermite_e
+
+This module provides a number of objects (mostly functions) useful for
+dealing with HermiteE series, including a `HermiteE` class that
+encapsulates the usual arithmetic operations.  (General information
+on how this module represents and works with such polynomials is in the
+docstring for its "parent" sub-package, `numpy.polynomial`).
+
+HermiteE Class
+--------------
+
+.. autosummary::
+   :toctree: generated/
+
+   HermiteE
+
+Basics
+------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermeval
+   hermeval2d
+   hermeval3d
+   hermegrid2d
+   hermegrid3d
+   hermeroots
+   hermefromroots
+
+Fitting
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermefit
+   hermevander
+   hermevander2d
+   hermevander3d
+
+Calculus
+--------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermeder
+   hermeint
+
+Algebra
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermeadd
+   hermesub
+   hermemul
+   hermemulx
+   hermediv
+   hermepow
+
+Quadrature
+----------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermegauss
+   hermeweight
+
+Miscellaneous
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   hermecompanion
+   hermedomain
+   hermezero
+   hermeone
+   hermex
+   hermetrim
+   hermeline
+   herme2poly
+   poly2herme
diff --git a/doc/source/reference/routines.polynomials.laguerre.rst b/doc/source/reference/routines.polynomials.laguerre.rst
new file mode 100644
index 000000000..45e288cb9
--- /dev/null
+++ b/doc/source/reference/routines.polynomials.laguerre.rst
@@ -0,0 +1,92 @@
+Laguerre Module (:mod:`numpy.polynomial.laguerre`)
+==================================================
+
+.. versionadded:: 1.6.0
+
+.. currentmodule:: numpy.polynomial.laguerre
+
+This module provides a number of objects (mostly functions) useful for
+dealing with Laguerre series, including a `Laguerre` class that
+encapsulates the usual arithmetic operations.  (General information
+on how this module represents and works with such polynomials is in the
+docstring for its "parent" sub-package, `numpy.polynomial`).
+
+Laguerre Class
+--------------
+
+.. autosummary::
+   :toctree: generated/
+
+   Laguerre
+
+Basics
+------
+
+.. autosummary::
+   :toctree: generated/
+
+   lagval
+   lagval2d
+   lagval3d
+   laggrid2d
+   laggrid3d
+   lagroots
+   lagfromroots
+
+Fitting
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   lagfit
+   lagvander
+   lagvander2d
+   lagvander3d
+
+Calculus
+--------
+
+.. autosummary::
+   :toctree: generated/
+
+   lagder
+   lagint
+
+Algebra
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   lagadd
+   lagsub
+   lagmul
+   lagmulx
+   lagdiv
+   lagpow
+
+Quadrature
+----------
+
+.. autosummary::
+   :toctree: generated/
+
+   laggauss
+   lagweight
+
+Miscellaneous
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   lagcompanion
+   lagdomain
+   lagzero
+   lagone
+   lagx
+   lagtrim
+   lagline
+   lag2poly
+   poly2lag
diff --git a/doc/source/reference/routines.polynomials.legendre.rst b/doc/source/reference/routines.polynomials.legendre.rst
new file mode 100644
index 000000000..fe6edc216
--- /dev/null
+++ b/doc/source/reference/routines.polynomials.legendre.rst
@@ -0,0 +1,92 @@
+Legendre Module (:mod:`numpy.polynomial.legendre`)
+==================================================
+
+.. versionadded:: 1.6.0
+
+.. currentmodule:: numpy.polynomial.legendre
+
+This module provides a number of objects (mostly functions) useful for
+dealing with Legendre series, including a `Legendre` class that
+encapsulates the usual arithmetic operations.  (General information
+on how this module represents and works with such polynomials is in the
+docstring for its "parent" sub-package, `numpy.polynomial`).
+
+Legendre Class
+--------------
+
+.. autosummary::
+   :toctree: generated/
+
+   Legendre
+
+Basics
+------
+
+.. autosummary::
+   :toctree: generated/
+
+   legval
+   legval2d
+   legval3d
+   leggrid2d
+   leggrid3d
+   legroots
+   legfromroots
+
+Fitting
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   legfit
+   legvander
+   legvander2d
+   legvander3d
+
+Calculus
+--------
+
+.. autosummary::
+   :toctree: generated/
+
+   legder
+   legint
+
+Algebra
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   legadd
+   legsub
+   legmul
+   legmulx
+   legdiv
+   legpow
+
+Quadrature
+----------
+
+.. autosummary::
+   :toctree: generated/
+
+   leggauss
+   legweight
+
+Miscellaneous
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   legcompanion
+   legdomain
+   legzero
+   legone
+   legx
+   legtrim
+   legline
+   leg2poly
+   poly2leg
diff --git a/doc/source/reference/routines.polynomials.package.rst b/doc/source/reference/routines.polynomials.package.rst
new file mode 100644
index 000000000..b2d357b31
--- /dev/null
+++ b/doc/source/reference/routines.polynomials.package.rst
@@ -0,0 +1,17 @@
+Polynomial Package
+==================
+
+.. versionadded:: 1.4.0
+
+.. currentmodule:: numpy.polynomial
+
+.. toctree::
+   :maxdepth: 2
+
+   routines.polynomials.classes
+   routines.polynomials.polynomial
+   routines.polynomials.chebyshev
+   routines.polynomials.legendre
+   routines.polynomials.laguerre
+   routines.polynomials.hermite
+   routines.polynomials.hermite_e
diff --git a/doc/source/reference/routines.polynomials.polynomial.rst b/doc/source/reference/routines.polynomials.polynomial.rst
index aa92ce8fc..431856622 100644
--- a/doc/source/reference/routines.polynomials.polynomial.rst
+++ b/doc/source/reference/routines.polynomials.polynomial.rst
@@ -1,16 +1,81 @@
-Polynomial Package (:mod:`numpy.polynomial`)
-============================================
+Polynomial Module (:mod:`numpy.polynomial.polynomial`)
+======================================================
 
-.. currentmodule:: numpy.polynomial
+.. versionadded:: 1.4.0
+
+.. currentmodule:: numpy.polynomial.polynomial
+
+This module provides a number of objects (mostly functions) useful for
+dealing with Polynomial series, including a `Polynomial` class that
+encapsulates the usual arithmetic operations.  (General information
+on how this module represents and works with such polynomials is in the
+docstring for its "parent" sub-package, `numpy.polynomial`).
+
+Polynomial Class
+----------------
 
-Polynomial Classes
-------------------
 .. autosummary::
    :toctree: generated/
 
    Polynomial
-   Chebyshev
-   Legendre
-   Hermite
-   HermiteE
-   Laguerre
+
+Basics
+------
+
+.. autosummary::
+   :toctree: generated/
+
+   polyval
+   polyval2d
+   polyval3d
+   polygrid2d
+   polygrid3d
+   polyroots
+   polyfromroots
+
+Fitting
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   polyfit
+   polyvander
+   polyvander2d
+   polyvander3d
+
+Calculus
+--------
+
+.. autosummary::
+   :toctree: generated/
+
+   polyder
+   polyint
+
+Algebra
+-------
+
+.. autosummary::
+   :toctree: generated/
+
+   polyadd
+   polysub
+   polymul
+   polymulx
+   polydiv
+   polypow
+
+Miscellaneous
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+   polycompanion
+   polydomain
+   polyzero
+   polyone
+   polyx
+   polytrim
+   polyline
diff --git a/doc/source/reference/routines.polynomials.rst b/doc/source/reference/routines.polynomials.rst
index 59d6bc499..94d1af8e7 100644
--- a/doc/source/reference/routines.polynomials.rst
+++ b/doc/source/reference/routines.polynomials.rst
@@ -1,14 +1,23 @@
 Polynomials
 ***********
 
-The poly1d functions are considered outdated but are retained for
-backward compatibility. New software needing polynomials should
-use the classes in the Polynomial Package.
+The polynomial package is newer and more complete than poly1d and the
+convenience classes are better behaved in the numpy environment. When
+backwards compatibility is not an issue it should be the package of choice.
+Note that the  various routines in the polynomial package all deal with
+series whose coefficients go from degree zero upward, which is the reverse
+of the poly1d convention. The easy way to remember this is that indexes
+correspond to degree, i.e., coef[i] is the coefficient of the term of
+degree i.
+
 
 .. toctree::
    :maxdepth: 2
 
-   routines.polynomials.polynomial
    routines.polynomials.poly1d
 
 
+.. toctree::
+   :maxdepth: 3
+
+   routines.polynomials.package
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 6b1edf497..e32cb229a 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -525,10 +525,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()
+        p = [chebline(-r, 1) for r in roots]
+        n = len(p)
+        while n > 1:
+            m, r = divmod(n, 2)
+            tmp = [chebmul(p[i], p[i+m]) for i in range(m)]
+            if r:
+                tmp[0] = chebmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def chebadd(c1, c2):
@@ -1122,9 +1129,16 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Chebyshev series to data.
 
-    Fit a Chebyshev series ``p(x) = p[0] * T_{0}(x) + ... + p[deg] *
-    T_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of
-    coefficients `p` that minimises the squared error.
+    Return the coefficients of a Legendre series of degree `deg` that is the
+    least squares fit to the data values `y` given at points `x`. If `y` is
+    1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
+    fits are done, one for each column of `y`, and the resulting
+    coefficients are stored in the corresponding columns of a 2-D return.
+    The fitted polynomial(s) are in the form
+
+    .. math::  p(x) = c_0 + c_1 * T_1(x) + ... + c_n * T_n(x),
+
+    where `n` is `deg`.
 
     Parameters
     ----------
@@ -1135,7 +1149,7 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
     deg : int
-        Degree of the fitting polynomial
+        Degree of the fitting series
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1150,6 +1164,7 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
         ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
         weights are chosen so that the errors of the products ``w[i]*y[i]``
         all have the same variance.  The default value is None.
+
         .. versionadded:: 1.5.0
 
     Returns
@@ -1176,30 +1191,31 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
 
     See Also
     --------
+    polyfit, legfit, lagfit, hermfit, hermefit
     chebval : Evaluates a Chebyshev series.
     chebvander : Vandermonde matrix of Chebyshev series.
-    polyfit : least squares fit using polynomials.
+    chebweight : Chebyshev weight function.
     linalg.lstsq : Computes a least-squares fit from the matrix.
     scipy.interpolate.UnivariateSpline : Computes spline fits.
 
     Notes
     -----
-    The solution are the coefficients ``c[i]`` of the Chebyshev series
-    ``T(x)`` that minimizes the squared error
+    The solution is the coefficients of the Chebyshev series `p` that
+    minimizes the sum of the weighted squared errors
 
-    ``E = \\sum_j |y_j - T(x_j)|^2``.
+    .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
 
-    This problem is solved by setting up as the overdetermined matrix
-    equation
+    where :math:`w_j` are the weights. This problem is solved by setting up
+    as the (typically) overdetermined matrix equation
 
-    ``V(x)*c = y``,
+    .. math:: V(x) * c = w * y,
 
-    where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are
-    the coefficients to be solved for, and the elements of `y` are the
+    where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
+    coefficients to be solved for, `w` are the weights, and `y` are the
     observed values.  This equation is then solved using the singular value
-    decomposition of ``V``.
+    decomposition of `V`.
 
-    If some of the singular values of ``V`` are so small that they are
+    If some of the singular values of `V` are so small that they are
     neglected, then a `RankWarning` will be issued. This means that the
     coeficient values may be poorly determined. Using a lower order fit
     will usually get rid of the warning.  The `rcond` parameter can also be
@@ -1273,6 +1289,46 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
         return c
 
 
+def chebcompanion(cs):
+    """Return the scaled companion matrix of cs.
+
+    The basis polynomials are scaled so that the companion matrix is
+    symmetric when `cs` represents a single Chebyshev polynomial. This
+    provides better eigenvalue estimates than the unscaled case and in the
+    single polynomial case the eigenvalues are guaranteed to be real if
+    np.eigvalsh is used to obtain them.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Legendre 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)
+    scl = np.array([1.] + [np.sqrt(.5)]*(n-1))
+    top = mat.reshape(-1)[1::n+1]
+    bot = mat.reshape(-1)[n::n+1]
+    top[0] = np.sqrt(.5)
+    top[1:] = 1/2
+    bot[...] = top
+    mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1])*.5
+    return mat
+
+
 def chebroots(cs):
     """
     Compute the roots of a Chebyshev series.
@@ -1307,34 +1363,22 @@ def chebroots(cs):
 
     Examples
     --------
-    >>> import numpy.polynomial as P
-    >>> import numpy.polynomial.chebyshev as C
-    >>> P.polyroots((-1,1,-1,1)) # x^3 - x^2 + x - 1 has two complex roots
-    array([ -4.99600361e-16-1.j,  -4.99600361e-16+1.j,   1.00000e+00+0.j])
-    >>> C.chebroots((-1,1,-1,1)) # T3 - T2 + T1 - T0 has only real roots
+    >>> import numpy.polynomial.chebyshev as cheb
+    >>> cheb.chebroots((-1, 1,-1, 1)) # T3 - T2 + T1 - T0 has real roots
     array([ -5.00000000e-01,   2.60860684e-17,   1.00000000e+00])
 
     """
     # 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
-    cs /= cs[-1]
-    cmat = np.zeros((n,n), dtype=cs.dtype)
-    cmat[1, 0] = 1
-    for i in range(1, n):
-        cmat[i - 1, i] = .5
-        if i != n - 1:
-            cmat[i + 1, i] = .5
-        else:
-            cmat[:, i] -= cs[:-1]*.5
-    roots = la.eigvals(cmat)
-    roots.sort()
-    return roots
+    m = chebcompanion(cs)
+    r = la.eigvals(m)
+    r.sort()
+    return r
 
 
 def chebpts1(npts):
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index d266a6453..199b47e15 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -277,10 +277,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()
+        p = [hermline(-r, 1) for r in roots]
+        n = len(p)
+        while n > 1:
+            m, r = divmod(n, 2)
+            tmp = [hermmul(p[i], p[i+m]) for i in range(m)]
+            if r:
+                tmp[0] = hermmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def hermadd(c1, c2):
@@ -919,9 +926,16 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Hermite series to data.
 
-    Fit a Hermite series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] *
-    P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of
-    coefficients `p` that minimises the squared error.
+    Return the coefficients of a Hermite series of degree `deg` that is the
+    least squares fit to the data values `y` given at points `x`. If `y` is
+    1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
+    fits are done, one for each column of `y`, and the resulting
+    coefficients are stored in the corresponding columns of a 2-D return.
+    The fitted polynomial(s) are in the form
+
+    .. math::  p(x) = c_0 + c_1 * H_1(x) + ... + c_n * H_n(x),
+
+    where `n` is `deg`.
 
     Parameters
     ----------
@@ -972,41 +986,42 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
 
     See Also
     --------
+    chebfit, legfit, lagfit, polyfit, hermefit
     hermval : Evaluates a Hermite series.
     hermvander : Vandermonde matrix of Hermite series.
-    polyfit : least squares fit using polynomials.
-    chebfit : least squares fit using Chebyshev series.
+    hermweight : Hermite weight function
     linalg.lstsq : Computes a least-squares fit from the matrix.
     scipy.interpolate.UnivariateSpline : Computes spline fits.
 
     Notes
     -----
-    The solution are the coefficients ``c[i]`` of the Hermite series
-    ``P(x)`` that minimizes the squared error
+    The solution is the coefficients of the Hermite series `p` that
+    minimizes the sum of the weighted squared errors
 
-    ``E = \\sum_j |y_j - P(x_j)|^2``.
+    .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
 
-    This problem is solved by setting up as the overdetermined matrix
-    equation
+    where the :math:`w_j` are the weights. This problem is solved by
+    setting up the (typically) overdetermined matrix equation
 
-    ``V(x)*c = y``,
+    .. math:: V(x) * c = w * y,
 
-    where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are
-    the coefficients to be solved for, and the elements of `y` are the
+    where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
+    coefficients to be solved for, `w` are the weights, `y` are the
     observed values.  This equation is then solved using the singular value
-    decomposition of ``V``.
+    decomposition of `V`.
 
-    If some of the singular values of ``V`` are so small that they are
+    If some of the singular values of `V` are so small that they are
     neglected, then a `RankWarning` will be issued. This means that the
     coeficient values may be poorly determined. Using a lower order fit
     will usually get rid of the warning.  The `rcond` parameter can also be
     set to a value smaller than its default, but the resulting fit may be
     spurious and have large contributions from roundoff error.
 
-    Fits using Hermite series are usually better conditioned than fits
-    using power series, but much can depend on the distribution of the
-    sample points and the smoothness of the data. If the quality of the fit
-    is inadequate splines may be a good alternative.
+    Fits using Hermite series are probably most useful when the data can be
+    approximated by ``sqrt(w(x)) * p(x)``, where `w(x)` is the Hermite
+    weight. In that case the wieght ``sqrt(w(x[i])`` should be used
+    together with data values ``y[i]/sqrt(w(x[i])``. The weight function is
+    available as `hermweight`.
 
     References
     ----------
@@ -1076,6 +1091,47 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
         return c
 
 
+def hermcompanion(cs):
+    """Return the scaled companion matrix of cs.
+
+    The basis polynomials are scaled so that the companion matrix is
+    symmetric when `cs` represents a single Hermite polynomial. This
+    provides better eigenvalue estimates than the unscaled case and in the
+    single polynomial case the eigenvalues are guaranteed to be real if
+    `numpy.linalg.eigvalsh` is used to obtain them.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Legendre series coefficients ordered from low to high
+        degree.
+
+    Returns
+    -------
+    mat : ndarray
+        Scaled companion matrix of dimensions (deg, deg).
+
+    """
+    accprod = np.multiply.accumulate
+    # 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(-.5*cs[0]/cs[1])
+
+    n = len(cs) - 1
+    mat = np.zeros((n, n), dtype=cs.dtype)
+    scl = np.hstack((1., np.sqrt(2.*np.arange(1,n))))
+    scl = np.multiply.accumulate(scl)
+    top = mat.reshape(-1)[1::n+1]
+    bot = mat.reshape(-1)[n::n+1]
+    top[...] = np.sqrt(.5*np.arange(1,n))
+    bot[...] = top
+    mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1])*.5
+    return mat
+
+
 def hermroots(cs):
     """
     Compute the roots of a Hermite series.
@@ -1125,19 +1181,10 @@ def hermroots(cs):
     if len(cs) == 2 :
         return np.array([-.5*cs[0]/cs[1]])
 
-    n = len(cs) - 1
-    cs /= cs[-1]
-    cmat = np.zeros((n,n), dtype=cs.dtype)
-    cmat[1, 0] = .5
-    for i in range(1, n):
-        cmat[i - 1, i] = i
-        if i != n - 1:
-            cmat[i + 1, i] = .5
-        else:
-            cmat[:, i] -= cs[:-1]*.5
-    roots = la.eigvals(cmat)
-    roots.sort()
-    return roots
+    m = hermcompanion(cs)
+    r = la.eigvals(m)
+    r.sort()
+    return r
 
 
 #
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index e644e345a..d0efe7209 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -277,10 +277,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()
+        p = [hermeline(-r, 1) for r in roots]
+        n = len(p)
+        while n > 1:
+            m, r = divmod(n, 2)
+            tmp = [hermemul(p[i], p[i+m]) for i in range(m)]
+            if r:
+                tmp[0] = hermemul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def hermeadd(c1, c2):
@@ -915,9 +922,16 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Hermite series to data.
 
-    Fit a Hermite series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] *
-    P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of
-    coefficients `p` that minimises the squared error.
+    Return the coefficients of a HermiteE series of degree `deg` that is
+    the least squares fit to the data values `y` given at points `x`. If
+    `y` is 1-D the returned coefficients will also be 1-D. If `y` is 2-D
+    multiple fits are done, one for each column of `y`, and the resulting
+    coefficients are stored in the corresponding columns of a 2-D return.
+    The fitted polynomial(s) are in the form
+
+    .. math::  p(x) = c_0 + c_1 * He_1(x) + ... + c_n * He_n(x),
+
+    where `n` is `deg`.
 
     Parameters
     ----------
@@ -968,41 +982,42 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
 
     See Also
     --------
+    chebfit, legfit, polyfit, hermfit, polyfit
     hermeval : Evaluates a Hermite series.
-    hermevander : Vandermonde matrix of Hermite series.
-    polyfit : least squares fit using polynomials.
-    chebfit : least squares fit using Chebyshev series.
+    hermevander : pseudo Vandermonde matrix of Hermite series.
+    hermeweight : HermiteE weight function.
     linalg.lstsq : Computes a least-squares fit from the matrix.
     scipy.interpolate.UnivariateSpline : Computes spline fits.
 
     Notes
     -----
-    The solution are the coefficients ``c[i]`` of the Hermite series
-    ``P(x)`` that minimizes the squared error
+    The solution is the coefficients of the HermiteE series `p` that
+    minimizes the sum of the weighted squared errors
 
-    ``E = \\sum_j |y_j - P(x_j)|^2``.
+    .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
 
-    This problem is solved by setting up as the overdetermined matrix
-    equation
+    where the :math:`w_j` are the weights. This problem is solved by
+    setting up the (typically) overdetermined matrix equation
 
-    ``V(x)*c = y``,
+    .. math:: V(x) * c = w * y,
 
-    where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are
-    the coefficients to be solved for, and the elements of `y` are the
+    where `V` is the pseudo Vandermonde matrix of `x`, the elements of `c`
+    are the coefficients to be solved for, and the elements of `y` are the
     observed values.  This equation is then solved using the singular value
-    decomposition of ``V``.
+    decomposition of `V`.
 
-    If some of the singular values of ``V`` are so small that they are
+    If some of the singular values of `V` are so small that they are
     neglected, then a `RankWarning` will be issued. This means that the
     coeficient values may be poorly determined. Using a lower order fit
     will usually get rid of the warning.  The `rcond` parameter can also be
     set to a value smaller than its default, but the resulting fit may be
     spurious and have large contributions from roundoff error.
 
-    Fits using Hermite series are usually better conditioned than fits
-    using power series, but much can depend on the distribution of the
-    sample points and the smoothness of the data. If the quality of the fit
-    is inadequate splines may be a good alternative.
+    Fits using HermiteE series are probably most useful when the data can
+    be approximated by ``sqrt(w(x)) * p(x)``, where `w(x)` is the HermiteE
+    weight. In that case the wieght ``sqrt(w(x[i])`` should be used
+    together with data values ``y[i]/sqrt(w(x[i])``. The weight function is
+    available as `hermeweight`.
 
     References
     ----------
@@ -1011,7 +1026,7 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
 
     Examples
     --------
-    >>> from numpy.polynomial.hermite_e import hermefit, hermeval
+    >>> from numpy.polynomial.hermite_e import hermefik, hermeval
     >>> x = np.linspace(-10, 10)
     >>> err = np.random.randn(len(x))/10
     >>> y = hermeval(x, [1, 2, 3]) + err
@@ -1072,18 +1087,59 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
         return c
 
 
+def hermecompanion(cs):
+    """Return the scaled companion matrix of cs.
+
+    The basis polynomials are scaled so that the companion matrix is
+    symmetric when `cs` represents a single HermiteE polynomial. This
+    provides better eigenvalue estimates than the unscaled case and in the
+    single polynomial case the eigenvalues are guaranteed to be real if
+    `numpy.linalg.eigvalsh` is used to obtain them.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Legendre series coefficients ordered from low to high
+        degree.
+
+    Returns
+    -------
+    mat : ndarray
+        Scaled companion matrix of dimensions (deg, deg).
+
+    """
+    accprod = np.multiply.accumulate
+    # 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)
+    scl = np.hstack((1., np.sqrt(np.arange(1,n))))
+    scl = np.multiply.accumulate(scl)
+    top = mat.reshape(-1)[1::n+1]
+    bot = mat.reshape(-1)[n::n+1]
+    top[...] = np.sqrt(np.arange(1,n))
+    bot[...] = top
+    mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1])
+    return mat
+
+
 def hermeroots(cs):
     """
     Compute the roots of a Hermite series.
 
-    Return the roots (a.k.a "zeros") of the Hermite series represented by
+    Return the roots (a.k.a "zeros") of the HermiteE series represented by
     `cs`, which is the sequence of coefficients from lowest order "term"
     to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``.
 
     Parameters
     ----------
     cs : array_like
-        1-d array of Hermite series coefficients ordered from low to high.
+        1-d array of HermiteE series coefficients ordered from low to high.
 
     Returns
     -------
@@ -1119,21 +1175,12 @@ def hermeroots(cs):
     if len(cs) <= 1 :
         return np.array([], dtype=cs.dtype)
     if len(cs) == 2 :
-        return np.array([-.5*cs[0]/cs[1]])
+        return np.array([-cs[0]/cs[1]])
 
-    n = len(cs) - 1
-    cs /= cs[-1]
-    cmat = np.zeros((n,n), dtype=cs.dtype)
-    cmat[1, 0] = 1
-    for i in range(1, n):
-        cmat[i - 1, i] = i
-        if i != n - 1:
-            cmat[i + 1, i] = 1
-        else:
-            cmat[:, i] -= cs[:-1]
-    roots = la.eigvals(cmat)
-    roots.sort()
-    return roots
+    m = hermecompanion(cs)
+    r = la.eigvals(m)
+    r.sort()
+    return r
 
 
 #
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index b6389bf63..7b518404f 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -274,10 +274,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()
+        p = [lagline(-r, 1) for r in roots]
+        n = len(p)
+        while n > 1:
+            m, r = divmod(n, 2)
+            tmp = [lagmul(p[i], p[i+m]) for i in range(m)]
+            if r:
+                tmp[0] = lagmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def lagadd(c1, c2):
@@ -916,9 +923,16 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Laguerre series to data.
 
-    Fit a Laguerre series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] *
-    P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of
-    coefficients `p` that minimises the squared error.
+    Return the coefficients of a Laguerre series of degree `deg` that is the
+    least squares fit to the data values `y` given at points `x`. If `y` is
+    1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
+    fits are done, one for each column of `y`, and the resulting
+    coefficients are stored in the corresponding columns of a 2-D return.
+    The fitted polynomial(s) are in the form
+
+    .. math::  p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x),
+
+    where `n` is `deg`.
 
     Parameters
     ----------
@@ -969,41 +983,42 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
 
     See Also
     --------
+    chebfit, legfit, polyfit, hermfit, hermefit
     lagval : Evaluates a Laguerre series.
-    lagvander : Vandermonde matrix of Laguerre series.
-    polyfit : least squares fit using polynomials.
-    chebfit : least squares fit using Chebyshev series.
+    lagvander : pseudo Vandermonde matrix of Laguerre series.
+    lagweight : Laguerre weight function.
     linalg.lstsq : Computes a least-squares fit from the matrix.
     scipy.interpolate.UnivariateSpline : Computes spline fits.
 
     Notes
     -----
-    The solution are the coefficients ``c[i]`` of the Laguerre series
-    ``P(x)`` that minimizes the squared error
+    The solution is the coefficients of the Laguerre series `p` that
+    minimizes the sum of the weighted squared errors
 
-    ``E = \\sum_j |y_j - P(x_j)|^2``.
+    .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
 
-    This problem is solved by setting up as the overdetermined matrix
-    equation
+    where the :math:`w_j` are the weights. This problem is solved by
+    setting up as the (typically) overdetermined matrix equation
 
-    ``V(x)*c = y``,
+    .. math:: V(x) * c = w * y,
 
-    where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are
-    the coefficients to be solved for, and the elements of `y` are the
+    where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
+    coefficients to be solved for, `w` are the weights, and `y` are the
     observed values.  This equation is then solved using the singular value
-    decomposition of ``V``.
+    decomposition of `V`.
 
-    If some of the singular values of ``V`` are so small that they are
+    If some of the singular values of `V` are so small that they are
     neglected, then a `RankWarning` will be issued. This means that the
     coeficient values may be poorly determined. Using a lower order fit
     will usually get rid of the warning.  The `rcond` parameter can also be
     set to a value smaller than its default, but the resulting fit may be
     spurious and have large contributions from roundoff error.
 
-    Fits using Laguerre series are usually better conditioned than fits
-    using power series, but much can depend on the distribution of the
-    sample points and the smoothness of the data. If the quality of the fit
-    is inadequate splines may be a good alternative.
+    Fits using Laguerre series are probably most useful when the data can
+    be approximated by ``sqrt(w(x)) * p(x)``, where `w(x)` is the Laguerre
+    weight. In that case the wieght ``sqrt(w(x[i])`` should be used
+    together with data values ``y[i]/sqrt(w(x[i])``. The weight function is
+    available as `lagweight`.
 
     References
     ----------
@@ -1073,6 +1088,45 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
         return c
 
 
+def lagcompanion(cs):
+    """Return the companion matrix of cs.
+
+    The unscaled companion matrix of the Laguerre polynomials is already
+    symmetric when `cs` represents a single Laguerre polynomial, so no
+    further scaling is needed.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Laguerre series coefficients ordered from low to high
+        degree.
+
+    Returns
+    -------
+    mat : ndarray
+        Companion matrix of dimensions (deg, deg).
+
+    """
+    accprod = np.multiply.accumulate
+    # 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(1 + cs[0]/cs[1])
+
+    n = len(cs) - 1
+    mat = np.zeros((n, n), dtype=cs.dtype)
+    top = mat.reshape(-1)[1::n+1]
+    mid = mat.reshape(-1)[0::n+1]
+    bot = mat.reshape(-1)[n::n+1]
+    top[...] = -np.arange(1,n)
+    mid[...] = 2.*np.arange(n) + 1.
+    bot[...] = top
+    mat[:,-1] += (cs[:-1]/cs[-1])*n
+    return mat
+
+
 def lagroots(cs):
     """
     Compute the roots of a Laguerre series.
@@ -1122,21 +1176,10 @@ def lagroots(cs):
     if len(cs) == 2 :
         return np.array([1 + cs[0]/cs[1]])
 
-    n = len(cs) - 1
-    cs /= cs[-1]
-    cmat = np.zeros((n,n), dtype=cs.dtype)
-    cmat[0, 0] = 1
-    cmat[1, 0] = -1
-    for i in range(1, n):
-        cmat[i - 1, i] = -i
-        cmat[i, i] = 2*i + 1
-        if i != n - 1:
-            cmat[i + 1, i] = -(i + 1)
-        else:
-            cmat[:, i] += cs[:-1]*(i + 1)
-    roots = la.eigvals(cmat)
-    roots.sort()
-    return roots
+    m = lagcompanion(cs)
+    r = la.eigvals(m)
+    r.sort()
+    return r
 
 
 #
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 3ea68a5e6..94098632c 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -1,5 +1,8 @@
 """
-Objects for dealing with Legendre series.
+Legendre Series (:mod: `numpy.polynomial.legendre`)
+===================================================
+
+.. currentmodule:: numpy.polynomial.polynomial
 
 This module provides a number of objects (mostly functions) useful for
 dealing with Legendre series, including a `Legendre` class that
@@ -9,44 +12,64 @@ docstring for its "parent" sub-package, `numpy.polynomial`).
 
 Constants
 ---------
-- `legdomain` -- Legendre series default domain, [-1,1].
-- `legzero` -- Legendre series that evaluates identically to 0.
-- `legone` -- Legendre series that evaluates identically to 1.
-- `legx` -- Legendre series for the identity map, ``f(x) = x``.
+
+.. autosummary::
+   :toctree: generated/
+
+   legdomain            Legendre series default domain, [-1,1].
+   legzero              Legendre series that evaluates identically to 0.
+   legone               Legendre series that evaluates identically to 1.
+   legx                 Legendre series for the identity map, ``f(x) = x``.
 
 Arithmetic
 ----------
-- `legmulx` -- multiply a Legendre series in ``P_i(x)`` by ``x``.
-- `legadd` -- add two Legendre series.
-- `legsub` -- subtract one Legendre series from another.
-- `legmul` -- multiply two Legendre series.
-- `legdiv` -- divide one Legendre series by another.
-- `legpow` -- raise a Legendre series to an positive integer power
-- `legval` -- evaluate a Legendre series at given points.
+
+.. autosummary::
+   :toctree: generated/
+
+   legmulx              multiply a Legendre series in P_i(x) by x.
+   legadd               add two Legendre series.
+   legsub               subtract one Legendre series from another.
+   legmul               multiply two Legendre series.
+   legdiv               divide one Legendre series by another.
+   legpow               raise a Legendre series to an positive integer power
+   legval               evaluate a Legendre series at given points.
 
 Calculus
 --------
-- `legder` -- differentiate a Legendre series.
-- `legint` -- integrate a Legendre series.
+
+.. autosummary::
+   :toctree: generated/
+
+   legder               differentiate a Legendre series.
+   legint               integrate a Legendre series.
 
 Misc Functions
 --------------
-- `legfromroots` -- create a Legendre series with specified roots.
-- `legroots` -- find the roots of a Legendre series.
-- `legvander` -- Vandermonde-like matrix for Legendre polynomials.
-- `legfit` -- least-squares fit returning a Legendre series.
-- `legtrim` -- trim leading coefficients from a Legendre series.
-- `legline` -- Legendre series representing given straight line.
-- `leg2poly` -- convert a Legendre series to a polynomial.
-- `poly2leg` -- convert a polynomial to a Legendre series.
+
+.. autosummary::
+   :toctree: generated/
+
+   legfromroots          create a Legendre series with specified roots.
+   legroots              find the roots of a Legendre series.
+   legvander             Vandermonde-like matrix for Legendre polynomials.
+   legfit                least-squares fit returning a Legendre series.
+   legtrim               trim leading coefficients from a Legendre series.
+   legline               Legendre series representing given straight line.
+   leg2poly              convert a Legendre series to a polynomial.
+   poly2leg              convert a polynomial to a Legendre series.
 
 Classes
 -------
-- `Legendre` -- A Legendre series class.
+    Legendre            A Legendre series class.
 
 See also
 --------
-`numpy.polynomial`
+numpy.polynomial.polynomial
+numpy.polynomial.chebyshev
+numpy.polynomial.laguerre
+numpy.polynomial.hermite
+numpy.polynomial.hermite_e
 
 """
 from __future__ import division
@@ -283,10 +306,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()
+        p = [legline(-r, 1) for r in roots]
+        n = len(p)
+        while n > 1:
+            m, r = divmod(n, 2)
+            tmp = [legmul(p[i], p[i+m]) for i in range(m)]
+            if r:
+                tmp[0] = legmul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def legadd(c1, c2):
@@ -919,9 +949,16 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least squares fit of Legendre series to data.
 
-    Fit a Legendre series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] *
-    P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of
-    coefficients `p` that minimises the squared error.
+    Return the coefficients of a Legendre series of degree `deg` that is the
+    least squares fit to the data values `y` given at points `x`. If `y` is
+    1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
+    fits are done, one for each column of `y`, and the resulting
+    coefficients are stored in the corresponding columns of a 2-D return.
+    The fitted polynomial(s) are in the form
+
+    .. math::  p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x),
+
+    where `n` is `deg`.
 
     Parameters
     ----------
@@ -948,6 +985,8 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
         weights are chosen so that the errors of the products ``w[i]*y[i]``
         all have the same variance.  The default value is None.
 
+        .. versionadded:: 1.5.0
+
     Returns
     -------
     coef : ndarray, shape (M,) or (M, K)
@@ -972,31 +1011,31 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
 
     See Also
     --------
+    chebfit, polyfit, lagfit, hermfit, hermefit
     legval : Evaluates a Legendre series.
     legvander : Vandermonde matrix of Legendre series.
-    polyfit : least squares fit using polynomials.
-    chebfit : least squares fit using Chebyshev series.
+    legweight : Legendre weight function (= 1).
     linalg.lstsq : Computes a least-squares fit from the matrix.
     scipy.interpolate.UnivariateSpline : Computes spline fits.
 
     Notes
     -----
-    The solution are the coefficients ``c[i]`` of the Legendre series
-    ``P(x)`` that minimizes the squared error
+    The solution is the coefficients of the Legendre series `p` that
+    minimizes the sum of the weighted squared errors
 
-    ``E = \\sum_j |y_j - P(x_j)|^2``.
+    .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
 
-    This problem is solved by setting up as the overdetermined matrix
-    equation
+    where :math:`w_j` are the weights. This problem is solved by setting up
+    as the (typically) overdetermined matrix equation
 
-    ``V(x)*c = y``,
+    .. math:: V(x) * c = w * y,
 
-    where ``V`` is the Vandermonde matrix of `x`, the elements of ``c`` are
-    the coefficients to be solved for, and the elements of `y` are the
+    where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
+    coefficients to be solved for, `w` are the weights, and `y` are the
     observed values.  This equation is then solved using the singular value
-    decomposition of ``V``.
+    decomposition of `V`.
 
-    If some of the singular values of ``V`` are so small that they are
+    If some of the singular values of `V` are so small that they are
     neglected, then a `RankWarning` will be issued. This means that the
     coeficient values may be poorly determined. Using a lower order fit
     will usually get rid of the warning.  The `rcond` parameter can also be
@@ -1070,11 +1109,50 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
         return c
 
 
+def legcompanion(cs):
+    """Return the scaled companion matrix of cs.
+
+    The basis polynomials are scaled so that the companion matrix is
+    symmetric when `cs` represents a single Legendre polynomial. This
+    provides better eigenvalue estimates than the unscaled case and in the
+    single polynomial case the eigenvalues are guaranteed to be real if
+    `numpy.linalg.eigvalsh` is used to obtain them.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Legendre 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)
+    scl = 1./np.sqrt(2*np.arange(n) + 1)
+    top = mat.reshape(-1)[1::n+1]
+    bot = mat.reshape(-1)[n::n+1]
+    top[...] = np.arange(1, n)*scl[:n-1]*scl[1:n]
+    bot[...] = top
+    mat[:,-1] -= (cs[:-1]/cs[-1])*(scl/scl[-1])*(n/(2*n - 1))
+    return mat
+
+
 def legroots(cs):
     """
     Compute the roots of a Legendre series.
 
-    Return the roots (a.k.a "zeros") of the Legendre series represented by
+    Returns the roots (a.k.a "zeros") of the Legendre series represented by
     `cs`, which is the sequence of coefficients from lowest order "term"
     to highest, e.g., [1,2,3] is the series ``L_0 + 2*L_1 + 3*L_2``.
 
@@ -1082,12 +1160,15 @@ def legroots(cs):
     ----------
     cs : array_like
         1-d array of Legendre series coefficients ordered from low to high.
+    maxiter : int, optional
+        Maximum number of iterations of Newton to use in refining the
+        roots.
 
     Returns
     -------
     out : ndarray
-        Array of the roots.  If all the roots are real, then so is the
-        dtype of ``out``; otherwise, ``out``'s dtype is complex.
+        Sorted array of the roots. If all the roots are real, then so is
+        the dtype of ``out``; otherwise, ``out``'s dtype is complex.
 
     See Also
     --------
@@ -1096,43 +1177,36 @@ def legroots(cs):
 
     Notes
     -----
-    Algorithm(s) used:
-
-    Remember: because the Legendre series basis set is different from the
-    "standard" basis set, the results of this function *may* not be what
-    one is expecting.
+    The root estimates are obtained as the eigenvalues of the companion
+    matrix, Roots far from the real interval [-1, 1] in the complex plane
+    may have large errors due to the numerical instability of the Lengendre
+    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
+    interval [-1, 1] can be improved by a few iterations of Newton's
+    method.
+
+    The Legendre series basis polynomials aren't powers of ``x`` so the
+    results of this function may seem unintuitive.
 
     Examples
     --------
-    >>> import numpy.polynomial as P
-    >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots
-    array([-0.60582959+0.j        , -0.07208521-0.63832674j,
-           -0.07208521+0.63832674j])
-    >>> P.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots
+    >>> import numpy.polynomial.legendre as leg
+    >>> leg.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots
     array([-0.85099543, -0.11407192,  0.51506735])
 
     """
     # 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
-    cs /= cs[-1]
-    cmat = np.zeros((n,n), dtype=cs.dtype)
-    cmat[1, 0] = 1
-    for i in range(1, n):
-        tmp = 2*i + 1
-        cmat[i - 1, i] = i/tmp
-        if i != n - 1:
-            cmat[i + 1, i] = (i + 1)/tmp
-        else:
-            cmat[:, i] -= cs[:-1]*(i + 1)/tmp
-    roots = la.eigvals(cmat)
-    roots.sort()
-    return roots
+    m = legcompanion(cs)
+    r = la.eigvals(m)
+    r.sort()
+    return r
 
 
 #
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 3efe25920..6f67814ad 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -169,10 +169,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()
+        p = [polyline(-r, 1) for r in roots]
+        n = len(p)
+        while n > 1:
+            m, r = divmod(n, 2)
+            tmp = [polymul(p[i], p[i+m]) for i in range(m)]
+            if r:
+                tmp[0] = polymul(tmp[0], p[-1])
+            p = tmp
+            n = m
+        return p[0]
 
 
 def polyadd(c1, c2):
@@ -696,10 +703,16 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     """
     Least-squares fit of a polynomial to data.
 
-    Fit a polynomial ``c0 + c1*x + c2*x**2 + ... + c[deg]*x**deg`` to
-    points (`x`, `y`).  Returns a 1-d (if `y` is 1-d) or 2-d (if `y` is 2-d)
-    array of coefficients representing, from lowest order term to highest,
-    the polynomial(s) which minimize the total square error.
+    Return the coefficients of a polynomial of degree `deg` that is the
+    least squares fit to the data values `y` given at points `x`. If `y` is
+    1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
+    fits are done, one for each column of `y`, and the resulting
+    coefficients are stored in the corresponding columns of a 2-D return.
+    The fitted polynomial(s) are in the form
+
+    .. math::  p(x) = c_0 + c_1 * x + ... + c_n * x^n,
+
+    where `n` is `deg`.
 
     Parameters
     ----------
@@ -708,7 +721,7 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     y : array_like, shape (`M`,) or (`M`, `K`)
         y-coordinates of the sample points.  Several sets of sample points
         sharing the same x-coordinates can be (independently) fit with one
-        call to `polyfit` by passing in for `y` a 2-d array that contains
+        call to `polyfit` by passing in for `y` a 2-D array that contains
         one data set per column.
     deg : int
         Degree of the polynomial(s) to be fit.
@@ -728,12 +741,13 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
         ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
         weights are chosen so that the errors of the products ``w[i]*y[i]``
         all have the same variance.  The default value is None.
+
         .. versionadded:: 1.5.0
 
     Returns
     -------
     coef : ndarray, shape (`deg` + 1,) or (`deg` + 1, `K`)
-        Polynomial coefficients ordered from low to high.  If `y` was 2-d,
+        Polynomial coefficients ordered from low to high.  If `y` was 2-D,
         the coefficients in column `k` of `coef` represent the polynomial
         fit to the data in `y`'s `k`-th column.
 
@@ -755,27 +769,27 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
 
     See Also
     --------
+    chebfit, legfit, lagfit, hermfit, hermefit
     polyval : Evaluates a polynomial.
     polyvander : Vandermonde matrix for powers.
-    chebfit : least squares fit using Chebyshev series.
     linalg.lstsq : Computes a least-squares fit from the matrix.
     scipy.interpolate.UnivariateSpline : Computes spline fits.
 
     Notes
     -----
-    The solutions are the coefficients ``c[i]`` of the polynomial ``P(x)``
-    that minimizes the total squared error:
+    The solution is the coefficients of the polynomial `p` that minimizes
+    the sum of the weighted squared errors
 
-    .. math :: E = \\sum_j (y_j - P(x_j))^2
+    .. math :: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
 
-    This problem is solved by setting up the (typically) over-determined
-    matrix equation:
+    where the :math:`w_j` are the weights. This problem is solved by
+    setting up the (typically) over-determined matrix equation:
 
-    .. math :: V(x)*c = y
+    .. math :: V(x) * c = w * y,
 
-    where `V` is the Vandermonde matrix of `x`, the elements of `c` are the
-    coefficients to be solved for, and the elements of `y` are the observed
-    values.  This equation is then solved using the singular value
+    where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
+    coefficients to be solved for, `w` are the weights, and `y` are the
+    observed values.  This equation is then solved using the singular value
     decomposition of `V`.
 
     If some of the singular values of `V` are so small that they are
@@ -790,10 +804,11 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     contributions from roundoff error.
 
     Polynomial fits using double precision tend to "fail" at about
-    (polynomial) degree 20. Fits using Chebyshev series are generally
-    better conditioned, but much can still depend on the distribution of
-    the sample points and the smoothness of the data.  If the quality of
-    the fit is inadequate, splines may be a good alternative.
+    (polynomial) degree 20. Fits using Chebyshev or Legendre series are
+    generally better conditioned, but much can still depend on the
+    distribution of the sample points and the smoothness of the data.  If
+    the quality of the fit is inadequate, splines may be a good
+    alternative.
 
     Examples
     --------
@@ -872,6 +887,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.
@@ -896,32 +941,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
 
 
 #
diff --git a/numpy/polynomial/polytemplate.py b/numpy/polynomial/polytemplate.py
index 9d5607467..e390d3d77 100644
--- a/numpy/polynomial/polytemplate.py
+++ b/numpy/polynomial/polytemplate.py
@@ -19,8 +19,9 @@ else:
 
 polytemplate = string.Template('''
 from __future__ import division
-REL_IMPORT polyutils as pu
 import numpy as np
+import warnings
+REL_IMPORT polyutils as pu
 
 class $name(pu.PolyBase) :
     """A $name series class.
@@ -147,21 +148,25 @@ class $name(pu.PolyBase) :
         """
         return np.all(self.window == other.window)
 
-    def has_samewindow(self, other):
-        """Check if windows match.
+    def has_sametype(self, other):
+        """Check if types match.
 
         Parameters
         ----------
-        other : class instance
-            The other class must have the ``window`` attribute.
+        other : object
+            Class instance.
 
         Returns
         -------
         bool : boolean
-            True if the windows are the same, False otherwise.
+            True if other is same class as self
+
+        Notes
+        -----
+        .. versionadded:: 1.7.0
 
         """
-        return np.all(self.window == other.window)
+        return isinstance(other, self.__class__)
 
     def __init__(self, coef, domain=$domain, window=$domain) :
         [coef, dom, win] = pu.as_series([coef, domain, window], trim=False)
@@ -220,11 +225,15 @@ class $name(pu.PolyBase) :
 
     def __add__(self, other) :
         """Returns sum"""
-        if isinstance(other, self.__class__) :
-            if self.has_samedomain(other) and self.has_samewindow(other):
+        if isinstance(other, pu.PolyBase):
+            if not self.has_sametype(other):
+                raise TypeError("Polynomial types differ")
+            elif not self.has_samedomain(other):
+                raise TypeError("Domains differ")
+            elif not self.has_samewindow(other):
+                raise TypeError("Windows differ")
+            else:
                 coef = ${nick}add(self.coef, other.coef)
-            else :
-                raise PolyDomainError()
         else :
             try :
                 coef = ${nick}add(self.coef, other)
@@ -234,11 +243,15 @@ class $name(pu.PolyBase) :
 
     def __sub__(self, other) :
         """Returns difference"""
-        if isinstance(other, self.__class__) :
-            if self.has_samedomain(other) and self.has_samewindow(other):
+        if isinstance(other, pu.PolyBase):
+            if not self.has_sametype(other):
+                raise TypeError("Polynomial types differ")
+            elif not self.has_samedomain(other):
+                raise TypeError("Domains differ")
+            elif not self.has_samewindow(other):
+                raise TypeError("Windows differ")
+            else:
                 coef = ${nick}sub(self.coef, other.coef)
-            else :
-                raise PolyDomainError()
         else :
             try :
                 coef = ${nick}sub(self.coef, other)
@@ -248,11 +261,15 @@ class $name(pu.PolyBase) :
 
     def __mul__(self, other) :
         """Returns product"""
-        if isinstance(other, self.__class__) :
-            if self.has_samedomain(other) and self.has_samewindow(other):
+        if isinstance(other, pu.PolyBase):
+            if not self.has_sametype(other):
+                raise TypeError("Polynomial types differ")
+            elif not self.has_samedomain(other):
+                raise TypeError("Domains differ")
+            elif not self.has_samewindow(other):
+                raise TypeError("Windows differ")
+            else:
                 coef = ${nick}mul(self.coef, other.coef)
-            else :
-                raise PolyDomainError()
         else :
             try :
                 coef = ${nick}mul(self.coef, other)
@@ -261,32 +278,31 @@ class $name(pu.PolyBase) :
         return self.__class__(coef, self.domain, self.window)
 
     def __div__(self, other):
-        # set to __floordiv__ /.
+        # set to __floordiv__,  /, for now.
         return self.__floordiv__(other)
 
     def __truediv__(self, other) :
         # there is no true divide if the rhs is not a scalar, although it
         # could return the first n elements of an infinite series.
         # It is hard to see where n would come from, though.
-        if isinstance(other, self.__class__) :
-            if len(other.coef) == 1 :
-                coef = div(self.coef, other.coef)
-            else :
-                return NotImplemented
-        elif np.isscalar(other) :
+        if np.isscalar(other) :
             # this might be overly restrictive
             coef = self.coef/other
+            return self.__class__(coef, self.domain, self.window)
         else :
             return NotImplemented
-        return self.__class__(coef, self.domain, self.window)
 
     def __floordiv__(self, other) :
         """Returns the quotient."""
-        if isinstance(other, self.__class__) :
-            if np.all(self.domain == other.domain) :
+        if isinstance(other, pu.PolyBase):
+            if not self.has_sametype(other):
+                raise TypeError("Polynomial types differ")
+            elif not self.has_samedomain(other):
+                raise TypeError("Domains differ")
+            elif not self.has_samewindow(other):
+                raise TypeError("Windows differ")
+            else:
                 quo, rem = ${nick}div(self.coef, other.coef)
-            else :
-                raise PolyDomainError()
         else :
             try :
                 quo, rem = ${nick}div(self.coef, other)
@@ -296,11 +312,15 @@ class $name(pu.PolyBase) :
 
     def __mod__(self, other) :
         """Returns the remainder."""
-        if isinstance(other, self.__class__) :
-            if self.has_samedomain(other) and self.has_samewindow(other):
+        if isinstance(other, pu.PolyBase):
+            if not self.has_sametype(other):
+                raise TypeError("Polynomial types differ")
+            elif not self.has_samedomain(other):
+                raise TypeError("Domains differ")
+            elif not self.has_samewindow(other):
+                raise TypeError("Windows differ")
+            else:
                 quo, rem = ${nick}div(self.coef, other.coef)
-            else :
-                raise PolyDomainError()
         else :
             try :
                 quo, rem = ${nick}div(self.coef, other)
@@ -311,10 +331,12 @@ class $name(pu.PolyBase) :
     def __divmod__(self, other) :
         """Returns quo, remainder"""
         if isinstance(other, self.__class__) :
-            if self.has_samedomain(other) and self.has_samewindow(other):
+            if not self.has_samedomain(other):
+                raise TypeError("Domains are not equal")
+            elif not self.has_samewindow(other):
+                raise TypeError("Windows are not equal")
+            else:
                 quo, rem = ${nick}div(self.coef, other.coef)
-            else :
-                raise PolyDomainError()
         else :
             try :
                 quo, rem = ${nick}div(self.coef, other)
diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py
index 6026e0e64..97f3a1c72 100644
--- a/numpy/polynomial/tests/test_hermite_e.py
+++ b/numpy/polynomial/tests/test_hermite_e.py
@@ -255,7 +255,7 @@ class TestMisc(TestCase) :
 
     def test_hermeroots(self) :
         assert_almost_equal(herme.hermeroots([1]), [])
-        assert_almost_equal(herme.hermeroots([1, 1]), [-.5])
+        assert_almost_equal(herme.hermeroots([1, 1]), [-1])
         for i in range(2,5) :
             tgt = np.linspace(-1, 1, i)
             res = herme.hermeroots(herme.hermefromroots(tgt))
diff --git a/numpy/polynomial/tests/test_printing.py b/numpy/polynomial/tests/test_printing.py
new file mode 100644
index 000000000..9803d931c
--- /dev/null
+++ b/numpy/polynomial/tests/test_printing.py
@@ -0,0 +1,81 @@
+import numpy.polynomial as poly
+from numpy.testing import TestCase, run_module_suite, assert_
+
+class test_str(TestCase):
+    def test_polynomial_str(self):
+        res = str(poly.Polynomial([0,1]))
+        tgt = 'poly([0., 1.])'
+        assert_(res, tgt)
+
+
+    def test_chebyshev_str(self):
+        res = str(poly.Chebyshev([0,1]))
+        tgt = 'leg([0., 1.])'
+        assert_(res, tgt)
+
+
+    def test_legendre_str(self):
+        res = str(poly.Legendre([0,1]))
+        tgt = 'leg([0., 1.])'
+        assert_(res, tgt)
+
+
+    def test_hermite_str(self):
+        res = str(poly.Hermite([0,1]))
+        tgt = 'herm([0., 1.])'
+        assert_(res, tgt)
+
+
+    def test_hermiteE_str(self):
+        res = str(poly.HermiteE([0,1]))
+        tgt = 'herme([0., 1.])'
+        assert_(res, tgt)
+
+
+    def test_laguerre_str(self):
+        res = str(poly.Laguerre([0,1]))
+        tgt = 'lag([0., 1.])'
+        assert_(res, tgt)
+
+
+class test_repr(TestCase):
+    def test_polynomial_str(self):
+        res = repr(poly.Polynomial([0,1]))
+        tgt = 'Polynomial([0., 1.])'
+        assert_(res, tgt)
+
+
+    def test_chebyshev_str(self):
+        res = repr(poly.Chebyshev([0,1]))
+        tgt = 'Chebyshev([0., 1.], [-1., 1.], [-1., 1.])'
+        assert_(res, tgt)
+
+
+    def test_legendre_repr(self):
+        res = repr(poly.Legendre([0,1]))
+        tgt = 'Legendre([0., 1.], [-1., 1.], [-1., 1.])'
+        assert_(res, tgt)
+
+
+    def test_hermite_repr(self):
+        res = repr(poly.Hermite([0,1]))
+        tgt = 'Hermite([0., 1.], [-1., 1.], [-1., 1.])'
+        assert_(res, tgt)
+
+
+    def test_hermiteE_repr(self):
+        res = repr(poly.HermiteE([0,1]))
+        tgt = 'HermiteE([0., 1.], [-1., 1.], [-1., 1.])'
+        assert_(res, tgt)
+
+
+    def test_laguerre_repr(self):
+        res = repr(poly.Laguerre([0,1]))
+        tgt = 'Laguerre([0., 1.], [0., 1.], [0., 1.])'
+        assert_(res, tgt)
+
+
+#
+
+if __name__ == "__main__":
+    run_module_suite()