Merge branch 'poly'

12 files changed, 5307 insertions, 178 deletions

diff --git a/numpy/polynomial/__init__.py b/numpy/polynomial/__init__.py
index 6b004acc4..32b39aa98 100644
--- a/numpy/polynomial/__init__.py
+++ b/numpy/polynomial/__init__.py
@@ -17,6 +17,9 @@ from polynomial import *
 from chebyshev import *
 from legendre import *
 from polyutils import *
+from hermite import Hermite
+from hermite_e import HermiteE
+from laguerre import Laguerre
 
 from numpy.testing import Tester
 test = Tester().test
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
new file mode 100644
index 000000000..d266a6453
--- /dev/null
+++ b/numpy/polynomial/hermite.py
@@ -0,0 +1,1147 @@
+"""
+Objects for dealing with Hermite series.
+
+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`).
+
+Constants
+---------
+- `hermdomain` -- Hermite series default domain, [-1,1].
+- `hermzero` -- Hermite series that evaluates identically to 0.
+- `hermone` -- Hermite series that evaluates identically to 1.
+- `hermx` -- Hermite series for the identity map, ``f(x) = x``.
+
+Arithmetic
+----------
+- `hermmulx` -- multiply a Hermite series in ``P_i(x)`` by ``x``.
+- `hermadd` -- add two Hermite series.
+- `hermsub` -- subtract one Hermite series from another.
+- `hermmul` -- multiply two Hermite series.
+- `hermdiv` -- divide one Hermite series by another.
+- `hermval` -- evaluate a Hermite series at given points.
+
+Calculus
+--------
+- `hermder` -- differentiate a Hermite series.
+- `hermint` -- integrate a Hermite series.
+
+Misc Functions
+--------------
+- `hermfromroots` -- create a Hermite series with specified roots.
+- `hermroots` -- find the roots of a Hermite series.
+- `hermvander` -- Vandermonde-like matrix for Hermite polynomials.
+- `hermfit` -- least-squares fit returning a Hermite series.
+- `hermtrim` -- trim leading coefficients from a Hermite series.
+- `hermline` -- Hermite series of given straight line.
+- `herm2poly` -- convert a Hermite series to a polynomial.
+- `poly2herm` -- convert a polynomial to a Hermite series.
+
+Classes
+-------
+- `Hermite` -- A Hermite series class.
+
+See also
+--------
+`numpy.polynomial`
+
+"""
+from __future__ import division
+
+__all__ = ['hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline',
+        'hermadd', 'hermsub', 'hermmulx', 'hermmul', 'hermdiv', 'hermval',
+        'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots',
+        'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite']
+
+import numpy as np
+import numpy.linalg as la
+import polyutils as pu
+import warnings
+from polytemplate import polytemplate
+
+hermtrim = pu.trimcoef
+
+def poly2herm(pol) :
+    """
+    poly2herm(pol)
+
+    Convert a polynomial to a Hermite series.
+
+    Convert an array representing the coefficients of a polynomial (relative
+    to the "standard" basis) ordered from lowest degree to highest, to an
+    array of the coefficients of the equivalent Hermite series, ordered
+    from lowest to highest degree.
+
+    Parameters
+    ----------
+    pol : array_like
+        1-d array containing the polynomial coefficients
+
+    Returns
+    -------
+    cs : ndarray
+        1-d array containing the coefficients of the equivalent Hermite
+        series.
+
+    See Also
+    --------
+    herm2poly
+
+    Notes
+    -----
+    The easy way to do conversions between polynomial basis sets
+    is to use the convert method of a class instance.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import poly2herme
+    >>> poly2herm(np.arange(4))
+    array([ 1.   ,  2.75 ,  0.5  ,  0.375])
+
+    """
+    [pol] = pu.as_series([pol])
+    deg = len(pol) - 1
+    res = 0
+    for i in range(deg, -1, -1) :
+        res = hermadd(hermmulx(res), pol[i])
+    return res
+
+
+def herm2poly(cs) :
+    """
+    Convert a Hermite series to a polynomial.
+
+    Convert an array representing the coefficients of a Hermite series,
+    ordered from lowest degree to highest, to an array of the coefficients
+    of the equivalent polynomial (relative to the "standard" basis) ordered
+    from lowest to highest degree.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array containing the Hermite series coefficients, ordered
+        from lowest order term to highest.
+
+    Returns
+    -------
+    pol : ndarray
+        1-d array containing the coefficients of the equivalent polynomial
+        (relative to the "standard" basis) ordered from lowest order term
+        to highest.
+
+    See Also
+    --------
+    poly2herm
+
+    Notes
+    -----
+    The easy way to do conversions between polynomial basis sets
+    is to use the convert method of a class instance.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import herm2poly
+    >>> herm2poly([ 1.   ,  2.75 ,  0.5  ,  0.375])
+    array([ 0.,  1.,  2.,  3.])
+
+    """
+    from polynomial import polyadd, polysub, polymulx
+
+    [cs] = pu.as_series([cs])
+    n = len(cs)
+    if n == 1:
+        return cs
+    if n == 2:
+        cs[1] *= 2
+        return cs
+    else:
+        c0 = cs[-2]
+        c1 = cs[-1]
+        # i is the current degree of c1
+        for i in range(n - 1, 1, -1) :
+            tmp = c0
+            c0 = polysub(cs[i - 2], c1*(2*(i - 1)))
+            c1 = polyadd(tmp, polymulx(c1)*2)
+        return polyadd(c0, polymulx(c1)*2)
+
+#
+# These are constant arrays are of integer type so as to be compatible
+# with the widest range of other types, such as Decimal.
+#
+
+# Hermite
+hermdomain = np.array([-1,1])
+
+# Hermite coefficients representing zero.
+hermzero = np.array([0])
+
+# Hermite coefficients representing one.
+hermone = np.array([1])
+
+# Hermite coefficients representing the identity x.
+hermx = np.array([0, 1/2])
+
+
+def hermline(off, scl) :
+    """
+    Hermite series whose graph is a straight line.
+
+
+
+    Parameters
+    ----------
+    off, scl : scalars
+        The specified line is given by ``off + scl*x``.
+
+    Returns
+    -------
+    y : ndarray
+        This module's representation of the Hermite series for
+        ``off + scl*x``.
+
+    See Also
+    --------
+    polyline, chebline
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermline, hermval
+    >>> hermval(0,hermline(3, 2))
+    3.0
+    >>> hermval(1,hermline(3, 2))
+    5.0
+
+    """
+    if scl != 0 :
+        return np.array([off,scl/2])
+    else :
+        return np.array([off])
+
+
+def hermfromroots(roots) :
+    """
+    Generate a Hermite series with the given roots.
+
+    Return the array of coefficients for the P-series whose roots (a.k.a.
+    "zeros") are given by *roots*.  The returned array of coefficients is
+    ordered from lowest order "term" to highest, and zeros of multiplicity
+    greater than one must be included in *roots* a number of times equal
+    to their multiplicity (e.g., if `2` is a root of multiplicity three,
+    then [2,2,2] must be in *roots*).
+
+    Parameters
+    ----------
+    roots : array_like
+        Sequence containing the roots.
+
+    Returns
+    -------
+    out : ndarray
+        1-d array of the Hermite series coefficients, ordered from low to
+        high.  If all roots are real, ``out.dtype`` is a float type;
+        otherwise, ``out.dtype`` is a complex type, even if all the
+        coefficients in the result are real (see Examples below).
+
+    See Also
+    --------
+    polyfromroots, chebfromroots
+
+    Notes
+    -----
+    What is returned are the :math:`c_i` such that:
+
+    .. math::
+
+        \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i])
+
+    where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Hermite
+    (basis) polynomial over the domain `[-1,1]`.  Note that, unlike
+    `polyfromroots`, due to the nature of the Hermite basis set, the
+    above identity *does not* imply :math:`c_n = 1` identically (see
+    Examples).
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermfromroots, hermval
+    >>> coef = hermfromroots((-1, 0, 1))
+    >>> hermval((-1, 0, 1), coef)
+    array([ 0.,  0.,  0.])
+    >>> coef = hermfromroots((-1j, 1j))
+    >>> hermval((-1j, 1j), coef)
+    array([ 0.+0.j,  0.+0.j])
+
+    """
+    if len(roots) == 0 :
+        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
+
+
+def hermadd(c1, c2):
+    """
+    Add one Hermite series to another.
+
+    Returns the sum of two Hermite series `c1` + `c2`.  The arguments
+    are sequences of coefficients ordered from lowest order term to
+    highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the Hermite series of their sum.
+
+    See Also
+    --------
+    hermsub, hermmul, hermdiv, hermpow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the sum of two Hermite series
+    is a Hermite series (without having to "reproject" the result onto
+    the basis set) so addition, just like that of "standard" polynomials,
+    is simply "component-wise."
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermadd
+    >>> hermadd([1, 2, 3], [1, 2, 3, 4])
+    array([ 2.,  4.,  6.,  4.])
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if len(c1) > len(c2) :
+        c1[:c2.size] += c2
+        ret = c1
+    else :
+        c2[:c1.size] += c1
+        ret = c2
+    return pu.trimseq(ret)
+
+
+def hermsub(c1, c2):
+    """
+    Subtract one Hermite series from another.
+
+    Returns the difference of two Hermite series `c1` - `c2`.  The
+    sequences of coefficients are from lowest order term to highest, i.e.,
+    [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Hermite series coefficients representing their difference.
+
+    See Also
+    --------
+    hermadd, hermmul, hermdiv, hermpow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the difference of two Hermite
+    series is a Hermite series (without having to "reproject" the result
+    onto the basis set) so subtraction, just like that of "standard"
+    polynomials, is simply "component-wise."
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermsub
+    >>> hermsub([1, 2, 3, 4], [1, 2, 3])
+    array([ 0.,  0.,  0.,  4.])
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if len(c1) > len(c2) :
+        c1[:c2.size] -= c2
+        ret = c1
+    else :
+        c2 = -c2
+        c2[:c1.size] += c1
+        ret = c2
+    return pu.trimseq(ret)
+
+
+def hermmulx(cs):
+    """Multiply a Hermite series by x.
+
+    Multiply the Hermite series `cs` by x, where x is the independent
+    variable.
+
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the result of the multiplication.
+
+    Notes
+    -----
+    The multiplication uses the recursion relationship for Hermite
+    polynomials in the form
+
+    .. math::
+
+    xP_i(x) = (P_{i + 1}(x)/2 + i*P_{i - 1}(x))
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermmulx
+    >>> hermmulx([1, 2, 3])
+    array([ 2. ,  6.5,  1. ,  1.5])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    # The zero series needs special treatment
+    if len(cs) == 1 and cs[0] == 0:
+        return cs
+
+    prd = np.empty(len(cs) + 1, dtype=cs.dtype)
+    prd[0] = cs[0]*0
+    prd[1] = cs[0]/2
+    for i in range(1, len(cs)):
+        prd[i + 1] = cs[i]/2
+        prd[i - 1] += cs[i]*i
+    return prd
+
+
+def hermmul(c1, c2):
+    """
+    Multiply one Hermite series by another.
+
+    Returns the product of two Hermite series `c1` * `c2`.  The arguments
+    are sequences of coefficients, from lowest order "term" to highest,
+    e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Hermite series coefficients representing their product.
+
+    See Also
+    --------
+    hermadd, hermsub, hermdiv, hermpow
+
+    Notes
+    -----
+    In general, the (polynomial) product of two C-series results in terms
+    that are not in the Hermite polynomial basis set.  Thus, to express
+    the product as a Hermite series, it is necessary to "re-project" the
+    product onto said basis set, which may produce "un-intuitive" (but
+    correct) results; see Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermmul
+    >>> hermmul([1, 2, 3], [0, 1, 2])
+    array([ 52.,  29.,  52.,   7.,   6.])
+
+    """
+    # s1, s2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+
+    if len(c1) > len(c2):
+        cs = c2
+        xs = c1
+    else:
+        cs = c1
+        xs = c2
+
+    if len(cs) == 1:
+        c0 = cs[0]*xs
+        c1 = 0
+    elif len(cs) == 2:
+        c0 = cs[0]*xs
+        c1 = cs[1]*xs
+    else :
+        nd = len(cs)
+        c0 = cs[-2]*xs
+        c1 = cs[-1]*xs
+        for i in range(3, len(cs) + 1) :
+            tmp = c0
+            nd =  nd - 1
+            c0 = hermsub(cs[-i]*xs, c1*(2*(nd - 1)))
+            c1 = hermadd(tmp, hermmulx(c1)*2)
+    return hermadd(c0, hermmulx(c1)*2)
+
+
+def hermdiv(c1, c2):
+    """
+    Divide one Hermite series by another.
+
+    Returns the quotient-with-remainder of two Hermite series
+    `c1` / `c2`.  The arguments are sequences of coefficients from lowest
+    order "term" to highest, e.g., [1,2,3] represents the series
+    ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    [quo, rem] : ndarrays
+        Of Hermite series coefficients representing the quotient and
+        remainder.
+
+    See Also
+    --------
+    hermadd, hermsub, hermmul, hermpow
+
+    Notes
+    -----
+    In general, the (polynomial) division of one Hermite series by another
+    results in quotient and remainder terms that are not in the Hermite
+    polynomial basis set.  Thus, to express these results as a Hermite
+    series, it is necessary to "re-project" the results onto the Hermite
+    basis set, which may produce "un-intuitive" (but correct) results; see
+    Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermdiv
+    >>> hermdiv([ 52.,  29.,  52.,   7.,   6.], [0, 1, 2])
+    (array([ 1.,  2.,  3.]), array([ 0.]))
+    >>> hermdiv([ 54.,  31.,  52.,   7.,   6.], [0, 1, 2])
+    (array([ 1.,  2.,  3.]), array([ 2.,  2.]))
+    >>> hermdiv([ 53.,  30.,  52.,   7.,   6.], [0, 1, 2])
+    (array([ 1.,  2.,  3.]), array([ 1.,  1.]))
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if c2[-1] == 0 :
+        raise ZeroDivisionError()
+
+    lc1 = len(c1)
+    lc2 = len(c2)
+    if lc1 < lc2 :
+        return c1[:1]*0, c1
+    elif lc2 == 1 :
+        return c1/c2[-1], c1[:1]*0
+    else :
+        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
+        rem = c1
+        for i in range(lc1 - lc2, - 1, -1):
+            p = hermmul([0]*i + [1], c2)
+            q = rem[-1]/p[-1]
+            rem = rem[:-1] - q*p[:-1]
+            quo[i] = q
+        return quo, pu.trimseq(rem)
+
+
+def hermpow(cs, pow, maxpower=16) :
+    """Raise a Hermite series to a power.
+
+    Returns the Hermite series `cs` raised to the power `pow`. The
+    arguement `cs` is a sequence of coefficients ordered from low to high.
+    i.e., [1,2,3] is the series  ``P_0 + 2*P_1 + 3*P_2.``
+
+    Parameters
+    ----------
+    cs : array_like
+        1d array of Hermite series coefficients ordered from low to
+        high.
+    pow : integer
+        Power to which the series will be raised
+    maxpower : integer, optional
+        Maximum power allowed. This is mainly to limit growth of the series
+        to umanageable size. Default is 16
+
+    Returns
+    -------
+    coef : ndarray
+        Hermite series of power.
+
+    See Also
+    --------
+    hermadd, hermsub, hermmul, hermdiv
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermpow
+    >>> hermpow([1, 2, 3], 2)
+    array([ 81.,  52.,  82.,  12.,   9.])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    power = int(pow)
+    if power != pow or power < 0 :
+        raise ValueError("Power must be a non-negative integer.")
+    elif maxpower is not None and power > maxpower :
+        raise ValueError("Power is too large")
+    elif power == 0 :
+        return np.array([1], dtype=cs.dtype)
+    elif power == 1 :
+        return cs
+    else :
+        # This can be made more efficient by using powers of two
+        # in the usual way.
+        prd = cs
+        for i in range(2, power + 1) :
+            prd = hermmul(prd, cs)
+        return prd
+
+
+def hermder(cs, m=1, scl=1) :
+    """
+    Differentiate a Hermite series.
+
+    Returns the series `cs` differentiated `m` times.  At each iteration the
+    result is multiplied by `scl` (the scaling factor is for use in a linear
+    change of variable).  The argument `cs` is the sequence of coefficients
+    from lowest order "term" to highest, e.g., [1,2,3] represents the series
+    ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    cs: array_like
+        1-d array of Hermite series coefficients ordered from low to high.
+    m : int, optional
+        Number of derivatives taken, must be non-negative. (Default: 1)
+    scl : scalar, optional
+        Each differentiation is multiplied by `scl`.  The end result is
+        multiplication by ``scl**m``.  This is for use in a linear change of
+        variable. (Default: 1)
+
+    Returns
+    -------
+    der : ndarray
+        Hermite series of the derivative.
+
+    See Also
+    --------
+    hermint
+
+    Notes
+    -----
+    In general, the result of differentiating a Hermite series does not
+    resemble the same operation on a power series. Thus the result of this
+    function may be "un-intuitive," albeit correct; see Examples section
+    below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermder
+    >>> hermder([ 1. ,  0.5,  0.5,  0.5])
+    array([ 1.,  2.,  3.])
+    >>> hermder([-0.5,  1./2.,  1./8.,  1./12.,  1./16.], m=2)
+    array([ 1.,  2.,  3.])
+
+    """
+    cnt = int(m)
+
+    if cnt != m:
+        raise ValueError, "The order of derivation must be integer"
+    if cnt < 0 :
+        raise ValueError, "The order of derivation must be non-negative"
+
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if cnt == 0:
+        return cs
+    elif cnt >= len(cs):
+        return cs[:1]*0
+    else :
+        for i in range(cnt):
+            n = len(cs) - 1
+            cs *= scl
+            der = np.empty(n, dtype=cs.dtype)
+            for j in range(n, 0, -1):
+                der[j - 1] = (2*j)*cs[j]
+            cs = der
+        return cs
+
+
+def hermint(cs, m=1, k=[], lbnd=0, scl=1):
+    """
+    Integrate a Hermite series.
+
+    Returns a Hermite series that is the Hermite series `cs`, integrated
+    `m` times from `lbnd` to `x`.  At each iteration the resulting series
+    is **multiplied** by `scl` and an integration constant, `k`, is added.
+    The scaling factor is for use in a linear change of variable.  ("Buyer
+    beware": note that, depending on what one is doing, one may want `scl`
+    to be the reciprocal of what one might expect; for more information,
+    see the Notes section below.)  The argument `cs` is a sequence of
+    coefficients, from lowest order Hermite series "term" to highest,
+    e.g., [1,2,3] represents the series :math:`P_0(x) + 2P_1(x) + 3P_2(x)`.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Hermite series coefficients, ordered from low to high.
+    m : int, optional
+        Order of integration, must be positive. (Default: 1)
+    k : {[], list, scalar}, optional
+        Integration constant(s).  The value of the first integral at
+        ``lbnd`` is the first value in the list, the value of the second
+        integral at ``lbnd`` is the second value, etc.  If ``k == []`` (the
+        default), all constants are set to zero.  If ``m == 1``, a single
+        scalar can be given instead of a list.
+    lbnd : scalar, optional
+        The lower bound of the integral. (Default: 0)
+    scl : scalar, optional
+        Following each integration the result is *multiplied* by `scl`
+        before the integration constant is added. (Default: 1)
+
+    Returns
+    -------
+    S : ndarray
+        Hermite series coefficients of the integral.
+
+    Raises
+    ------
+    ValueError
+        If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or
+        ``np.isscalar(scl) == False``.
+
+    See Also
+    --------
+    hermder
+
+    Notes
+    -----
+    Note that the result of each integration is *multiplied* by `scl`.
+    Why is this important to note?  Say one is making a linear change of
+    variable :math:`u = ax + b` in an integral relative to `x`.  Then
+    :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a`
+    - perhaps not what one would have first thought.
+
+    Also note that, in general, the result of integrating a C-series needs
+    to be "re-projected" onto the C-series basis set.  Thus, typically,
+    the result of this function is "un-intuitive," albeit correct; see
+    Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermint
+    >>> hermint([1,2,3]) # integrate once, value 0 at 0.
+    array([ 1. ,  0.5,  0.5,  0.5])
+    >>> hermint([1,2,3], m=2) # integrate twice, value & deriv 0 at 0 
+    array([-0.5       ,  0.5       ,  0.125     ,  0.08333333,  0.0625    ])
+    >>> hermint([1,2,3], k=1) # integrate once, value 1 at 0.
+    array([ 2. ,  0.5,  0.5,  0.5])
+    >>> hermint([1,2,3], lbnd=-1) # integrate once, value 0 at -1
+    array([-2. ,  0.5,  0.5,  0.5])
+    >>> hermint([1,2,3], m=2, k=[1,2], lbnd=-1)
+    array([ 1.66666667, -0.5       ,  0.125     ,  0.08333333,  0.0625    ])
+
+    """
+    cnt = int(m)
+    if np.isscalar(k) :
+        k = [k]
+
+    if cnt != m:
+        raise ValueError, "The order of integration must be integer"
+    if cnt < 0 :
+        raise ValueError, "The order of integration must be non-negative"
+    if len(k) > cnt :
+        raise ValueError, "Too many integration constants"
+
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if cnt == 0:
+        return cs
+
+    k = list(k) + [0]*(cnt - len(k))
+    for i in range(cnt) :
+        n = len(cs)
+        cs *= scl
+        if n == 1 and cs[0] == 0:
+            cs[0] += k[i]
+        else:
+            tmp = np.empty(n + 1, dtype=cs.dtype)
+            tmp[0] = cs[0]*0
+            tmp[1] = cs[0]/2
+            for j in range(1, n):
+                tmp[j + 1] = cs[j]/(2*(j + 1))
+            tmp[0] += k[i] - hermval(lbnd, tmp)
+            cs = tmp
+    return cs
+
+
+def hermval(x, cs):
+    """Evaluate a Hermite series.
+
+    If `cs` is of length `n`, this function returns :
+
+    ``p(x) = cs[0]*P_0(x) + cs[1]*P_1(x) + ... + cs[n-1]*P_{n-1}(x)``
+
+    If x is a sequence or array then p(x) will have the same shape as x.
+    If r is a ring_like object that supports multiplication and addition
+    by the values in `cs`, then an object of the same type is returned.
+
+    Parameters
+    ----------
+    x : array_like, ring_like
+        Array of numbers or objects that support multiplication and
+        addition with themselves and with the elements of `cs`.
+    cs : array_like
+        1-d array of Hermite coefficients ordered from low to high.
+
+    Returns
+    -------
+    values : ndarray, ring_like
+        If the return is an ndarray then it has the same shape as `x`.
+
+    See Also
+    --------
+    hermfit
+
+    Examples
+    --------
+
+    Notes
+    -----
+    The evaluation uses Clenshaw recursion, aka synthetic division.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermval
+    >>> coef = [1,2,3]
+    >>> hermval(1, coef)
+    11.0
+    >>> hermval([[1,2],[3,4]], coef)
+    array([[  11.,   51.],
+           [ 115.,  203.]])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if isinstance(x, tuple) or isinstance(x, list) :
+        x = np.asarray(x)
+
+    x2 = x*2
+    if len(cs) == 1 :
+        c0 = cs[0]
+        c1 = 0
+    elif len(cs) == 2 :
+        c0 = cs[0]
+        c1 = cs[1]
+    else :
+        nd = len(cs)
+        c0 = cs[-2]
+        c1 = cs[-1]
+        for i in range(3, len(cs) + 1) :
+            tmp = c0
+            nd =  nd - 1
+            c0 = cs[-i] - c1*(2*(nd - 1))
+            c1 = tmp + c1*x2
+    return c0 + c1*x2
+
+
+def hermvander(x, deg) :
+    """Vandermonde matrix of given degree.
+
+    Returns the Vandermonde matrix of degree `deg` and sample points `x`.
+    This isn't a true Vandermonde matrix because `x` can be an arbitrary
+    ndarray and the Hermite polynomials aren't powers. If ``V`` is the
+    returned matrix and `x` is a 2d array, then the elements of ``V`` are
+    ``V[i,j,k] = P_k(x[i,j])``, where ``P_k`` is the Hermite polynomial
+    of degree ``k``.
+
+    Parameters
+    ----------
+    x : array_like
+        Array of points. The values are converted to double or complex
+        doubles. If x is scalar it is converted to a 1D array.
+    deg : integer
+        Degree of the resulting matrix.
+
+    Returns
+    -------
+    vander : Vandermonde matrix.
+        The shape of the returned matrix is ``x.shape + (deg+1,)``. The last
+        index is the degree.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermvander
+    >>> x = np.array([-1, 0, 1])
+    >>> hermvander(x, 3)
+    array([[ 1., -2.,  2.,  4.],
+           [ 1.,  0., -2., -0.],
+           [ 1.,  2.,  2., -4.]])
+
+    """
+    ideg = int(deg)
+    if ideg != deg:
+        raise ValueError("deg must be integer")
+    if ideg < 0:
+        raise ValueError("deg must be non-negative")
+
+    x = np.array(x, copy=0, ndmin=1) + 0.0
+    v = np.empty((ideg + 1,) + x.shape, dtype=x.dtype)
+    v[0] = x*0 + 1
+    if ideg > 0 :
+        x2 = x*2
+        v[1] = x2
+        for i in range(2, ideg + 1) :
+            v[i] = (v[i-1]*x2 - v[i-2]*(2*(i - 1)))
+    return np.rollaxis(v, 0, v.ndim)
+
+
+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.
+
+    Parameters
+    ----------
+    x : array_like, shape (M,)
+        x-coordinates of the M sample points ``(x[i], y[i])``.
+    y : array_like, shape (M,) or (M, K)
+        y-coordinates of the sample points. Several data sets of sample
+        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
+    rcond : float, optional
+        Relative condition number of the fit. Singular values smaller than
+        this relative to the largest singular value will be ignored. The
+        default value is len(x)*eps, where eps is the relative precision of
+        the float type, about 2e-16 in most cases.
+    full : bool, optional
+        Switch determining nature of return value. When it is False (the
+        default) just the coefficients are returned, when True diagnostic
+        information from the singular value decomposition is also returned.
+    w : array_like, shape (`M`,), optional
+        Weights. If not None, the contribution of each point
+        ``(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.
+
+    Returns
+    -------
+    coef : ndarray, shape (M,) or (M, K)
+        Hermite coefficients ordered from low to high. If `y` was 2-D,
+        the coefficients for the data in column k  of `y` are in column
+        `k`.
+
+    [residuals, rank, singular_values, rcond] : present when `full` = True
+        Residuals of the least-squares fit, the effective rank of the
+        scaled Vandermonde matrix and its singular values, and the
+        specified value of `rcond`. For more details, see `linalg.lstsq`.
+
+    Warns
+    -----
+    RankWarning
+        The rank of the coefficient matrix in the least-squares fit is
+        deficient. The warning is only raised if `full` = False.  The
+        warnings can be turned off by
+
+        >>> import warnings
+        >>> warnings.simplefilter('ignore', RankWarning)
+
+    See Also
+    --------
+    hermval : Evaluates a Hermite series.
+    hermvander : Vandermonde matrix of Hermite series.
+    polyfit : least squares fit using polynomials.
+    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 solution are the coefficients ``c[i]`` of the Hermite series
+    ``P(x)`` that minimizes the squared error
+
+    ``E = \\sum_j |y_j - P(x_j)|^2``.
+
+    This problem is solved by setting up as the overdetermined matrix
+    equation
+
+    ``V(x)*c = 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
+    decomposition of ``V``.
+
+    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.
+
+    References
+    ----------
+    .. [1] Wikipedia, "Curve fitting",
+           http://en.wikipedia.org/wiki/Curve_fitting
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermfit, hermval
+    >>> x = np.linspace(-10, 10)
+    >>> err = np.random.randn(len(x))/10
+    >>> y = hermval(x, [1, 2, 3]) + err
+    >>> hermfit(x, y, 2)
+    array([ 0.97902637,  1.99849131,  3.00006   ])
+
+    """
+    order = int(deg) + 1
+    x = np.asarray(x) + 0.0
+    y = np.asarray(y) + 0.0
+
+    # check arguments.
+    if deg < 0 :
+        raise ValueError, "expected deg >= 0"
+    if x.ndim != 1:
+        raise TypeError, "expected 1D vector for x"
+    if x.size == 0:
+        raise TypeError, "expected non-empty vector for x"
+    if y.ndim < 1 or y.ndim > 2 :
+        raise TypeError, "expected 1D or 2D array for y"
+    if len(x) != len(y):
+        raise TypeError, "expected x and y to have same length"
+
+    # set up the least squares matrices
+    lhs = hermvander(x, deg)
+    rhs = y
+    if w is not None:
+        w = np.asarray(w) + 0.0
+        if w.ndim != 1:
+            raise TypeError, "expected 1D vector for w"
+        if len(x) != len(w):
+            raise TypeError, "expected x and w to have same length"
+        # apply weights
+        if rhs.ndim == 2:
+            lhs *= w[:, np.newaxis]
+            rhs *= w[:, np.newaxis]
+        else:
+            lhs *= w[:, np.newaxis]
+            rhs *= w
+
+    # set rcond
+    if rcond is None :
+        rcond = len(x)*np.finfo(x.dtype).eps
+
+    # scale the design matrix and solve the least squares equation
+    scl = np.sqrt((lhs*lhs).sum(0))
+    c, resids, rank, s = la.lstsq(lhs/scl, rhs, rcond)
+    c = (c.T/scl).T
+
+    # warn on rank reduction
+    if rank != order and not full:
+        msg = "The fit may be poorly conditioned"
+        warnings.warn(msg, pu.RankWarning)
+
+    if full :
+        return c, [resids, rank, s, rcond]
+    else :
+        return c
+
+
+def hermroots(cs):
+    """
+    Compute the roots of a Hermite series.
+
+    Return the roots (a.k.a "zeros") of the Hermite 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.
+
+    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.
+
+    See Also
+    --------
+    polyroots
+    chebroots
+
+    Notes
+    -----
+    Algorithm(s) used:
+
+    Remember: because the Hermite series basis set is different from the
+    "standard" basis set, the results of this function *may* not be what
+    one is expecting.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite import hermroots, hermfromroots
+    >>> coef = hermfromroots([-1, 0, 1])
+    >>> coef
+    array([ 0.   ,  0.25 ,  0.   ,  0.125])
+    >>> hermroots(coef)
+    array([ -1.00000000e+00,  -1.38777878e-17,   1.00000000e+00])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) <= 1 :
+        return np.array([], dtype=cs.dtype)
+    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
+
+
+#
+# Hermite series class
+#
+
+exec polytemplate.substitute(name='Hermite', nick='herm', domain='[-1,1]')
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
new file mode 100644
index 000000000..e644e345a
--- /dev/null
+++ b/numpy/polynomial/hermite_e.py
@@ -0,0 +1,1143 @@
+"""
+Objects for dealing with Hermite series.
+
+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`).
+
+Constants
+---------
+- `hermedomain` -- Hermite series default domain, [-1,1].
+- `hermezero` -- Hermite series that evaluates identically to 0.
+- `hermeone` -- Hermite series that evaluates identically to 1.
+- `hermex` -- Hermite series for the identity map, ``f(x) = x``.
+
+Arithmetic
+----------
+- `hermemulx` -- multiply a Hermite series in ``P_i(x)`` by ``x``.
+- `hermeadd` -- add two Hermite series.
+- `hermesub` -- subtract one Hermite series from another.
+- `hermemul` -- multiply two Hermite series.
+- `hermediv` -- divide one Hermite series by another.
+- `hermeval` -- evaluate a Hermite series at given points.
+
+Calculus
+--------
+- `hermeder` -- differentiate a Hermite series.
+- `hermeint` -- integrate a Hermite series.
+
+Misc Functions
+--------------
+- `hermefromroots` -- create a Hermite series with specified roots.
+- `hermeroots` -- find the roots of a Hermite series.
+- `hermevander` -- Vandermonde-like matrix for Hermite polynomials.
+- `hermefit` -- least-squares fit returning a Hermite series.
+- `hermetrim` -- trim leading coefficients from a Hermite series.
+- `hermeline` -- Hermite series of given straight line.
+- `herme2poly` -- convert a Hermite series to a polynomial.
+- `poly2herme` -- convert a polynomial to a Hermite series.
+
+Classes
+-------
+- `Hermite` -- A Hermite series class.
+
+See also
+--------
+`numpy.polynomial`
+
+"""
+from __future__ import division
+
+import numpy as np
+import numpy.linalg as la
+import polyutils as pu
+import warnings
+from polytemplate import polytemplate
+
+__all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
+        'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermeval',
+        'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots',
+        'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE']
+
+hermetrim = pu.trimcoef
+
+def poly2herme(pol) :
+    """
+    poly2herme(pol)
+
+    Convert a polynomial to a Hermite series.
+
+    Convert an array representing the coefficients of a polynomial (relative
+    to the "standard" basis) ordered from lowest degree to highest, to an
+    array of the coefficients of the equivalent Hermite series, ordered
+    from lowest to highest degree.
+
+    Parameters
+    ----------
+    pol : array_like
+        1-d array containing the polynomial coefficients
+
+    Returns
+    -------
+    cs : ndarray
+        1-d array containing the coefficients of the equivalent Hermite
+        series.
+
+    See Also
+    --------
+    herme2poly
+
+    Notes
+    -----
+    The easy way to do conversions between polynomial basis sets
+    is to use the convert method of a class instance.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import poly2herme
+    >>> poly2herme(np.arange(4))
+    array([  2.,  10.,   2.,   3.])
+
+    """
+    [pol] = pu.as_series([pol])
+    deg = len(pol) - 1
+    res = 0
+    for i in range(deg, -1, -1) :
+        res = hermeadd(hermemulx(res), pol[i])
+    return res
+
+
+def herme2poly(cs) :
+    """
+    Convert a Hermite series to a polynomial.
+
+    Convert an array representing the coefficients of a Hermite series,
+    ordered from lowest degree to highest, to an array of the coefficients
+    of the equivalent polynomial (relative to the "standard" basis) ordered
+    from lowest to highest degree.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array containing the Hermite series coefficients, ordered
+        from lowest order term to highest.
+
+    Returns
+    -------
+    pol : ndarray
+        1-d array containing the coefficients of the equivalent polynomial
+        (relative to the "standard" basis) ordered from lowest order term
+        to highest.
+
+    See Also
+    --------
+    poly2herme
+
+    Notes
+    -----
+    The easy way to do conversions between polynomial basis sets
+    is to use the convert method of a class instance.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import herme2poly
+    >>> herme2poly([  2.,  10.,   2.,   3.])
+    array([ 0.,  1.,  2.,  3.])
+
+    """
+    from polynomial import polyadd, polysub, polymulx
+
+    [cs] = pu.as_series([cs])
+    n = len(cs)
+    if n == 1:
+        return cs
+    if n == 2:
+        return cs
+    else:
+        c0 = cs[-2]
+        c1 = cs[-1]
+        # i is the current degree of c1
+        for i in range(n - 1, 1, -1) :
+            tmp = c0
+            c0 = polysub(cs[i - 2], c1*(i - 1))
+            c1 = polyadd(tmp, polymulx(c1))
+        return polyadd(c0, polymulx(c1))
+
+#
+# These are constant arrays are of integer type so as to be compatible
+# with the widest range of other types, such as Decimal.
+#
+
+# Hermite
+hermedomain = np.array([-1,1])
+
+# Hermite coefficients representing zero.
+hermezero = np.array([0])
+
+# Hermite coefficients representing one.
+hermeone = np.array([1])
+
+# Hermite coefficients representing the identity x.
+hermex = np.array([0, 1])
+
+
+def hermeline(off, scl) :
+    """
+    Hermite series whose graph is a straight line.
+
+
+
+    Parameters
+    ----------
+    off, scl : scalars
+        The specified line is given by ``off + scl*x``.
+
+    Returns
+    -------
+    y : ndarray
+        This module's representation of the Hermite series for
+        ``off + scl*x``.
+
+    See Also
+    --------
+    polyline, chebline
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermeline
+    >>> from numpy.polynomial.hermite_e import hermeline, hermeval
+    >>> hermeval(0,hermeline(3, 2))
+    3.0
+    >>> hermeval(1,hermeline(3, 2))
+    5.0
+
+    """
+    if scl != 0 :
+        return np.array([off,scl])
+    else :
+        return np.array([off])
+
+
+def hermefromroots(roots) :
+    """
+    Generate a Hermite series with the given roots.
+
+    Return the array of coefficients for the P-series whose roots (a.k.a.
+    "zeros") are given by *roots*.  The returned array of coefficients is
+    ordered from lowest order "term" to highest, and zeros of multiplicity
+    greater than one must be included in *roots* a number of times equal
+    to their multiplicity (e.g., if `2` is a root of multiplicity three,
+    then [2,2,2] must be in *roots*).
+
+    Parameters
+    ----------
+    roots : array_like
+        Sequence containing the roots.
+
+    Returns
+    -------
+    out : ndarray
+        1-d array of the Hermite series coefficients, ordered from low to
+        high.  If all roots are real, ``out.dtype`` is a float type;
+        otherwise, ``out.dtype`` is a complex type, even if all the
+        coefficients in the result are real (see Examples below).
+
+    See Also
+    --------
+    polyfromroots, chebfromroots
+
+    Notes
+    -----
+    What is returned are the :math:`c_i` such that:
+
+    .. math::
+
+        \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i])
+
+    where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Hermite
+    (basis) polynomial over the domain `[-1,1]`.  Note that, unlike
+    `polyfromroots`, due to the nature of the Hermite basis set, the
+    above identity *does not* imply :math:`c_n = 1` identically (see
+    Examples).
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermefromroots, hermeval
+    >>> coef = hermefromroots((-1, 0, 1))
+    >>> hermeval((-1, 0, 1), coef)
+    array([ 0.,  0.,  0.])
+    >>> coef = hermefromroots((-1j, 1j))
+    >>> hermeval((-1j, 1j), coef)
+    array([ 0.+0.j,  0.+0.j])
+
+    """
+    if len(roots) == 0 :
+        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
+
+
+def hermeadd(c1, c2):
+    """
+    Add one Hermite series to another.
+
+    Returns the sum of two Hermite series `c1` + `c2`.  The arguments
+    are sequences of coefficients ordered from lowest order term to
+    highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the Hermite series of their sum.
+
+    See Also
+    --------
+    hermesub, hermemul, hermediv, hermepow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the sum of two Hermite series
+    is a Hermite series (without having to "reproject" the result onto
+    the basis set) so addition, just like that of "standard" polynomials,
+    is simply "component-wise."
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermeadd
+    >>> hermeadd([1, 2, 3], [1, 2, 3, 4])
+    array([ 2.,  4.,  6.,  4.])
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if len(c1) > len(c2) :
+        c1[:c2.size] += c2
+        ret = c1
+    else :
+        c2[:c1.size] += c1
+        ret = c2
+    return pu.trimseq(ret)
+
+
+def hermesub(c1, c2):
+    """
+    Subtract one Hermite series from another.
+
+    Returns the difference of two Hermite series `c1` - `c2`.  The
+    sequences of coefficients are from lowest order term to highest, i.e.,
+    [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Hermite series coefficients representing their difference.
+
+    See Also
+    --------
+    hermeadd, hermemul, hermediv, hermepow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the difference of two Hermite
+    series is a Hermite series (without having to "reproject" the result
+    onto the basis set) so subtraction, just like that of "standard"
+    polynomials, is simply "component-wise."
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermesub
+    >>> hermesub([1, 2, 3, 4], [1, 2, 3])
+    array([ 0.,  0.,  0.,  4.])
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if len(c1) > len(c2) :
+        c1[:c2.size] -= c2
+        ret = c1
+    else :
+        c2 = -c2
+        c2[:c1.size] += c1
+        ret = c2
+    return pu.trimseq(ret)
+
+
+def hermemulx(cs):
+    """Multiply a Hermite series by x.
+
+    Multiply the Hermite series `cs` by x, where x is the independent
+    variable.
+
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the result of the multiplication.
+
+    Notes
+    -----
+    The multiplication uses the recursion relationship for Hermite
+    polynomials in the form
+
+    .. math::
+
+    xP_i(x) = (P_{i + 1}(x) + iP_{i - 1}(x)))
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermemulx
+    >>> hermemulx([1, 2, 3])
+    array([ 2.,  7.,  2.,  3.])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    # The zero series needs special treatment
+    if len(cs) == 1 and cs[0] == 0:
+        return cs
+
+    prd = np.empty(len(cs) + 1, dtype=cs.dtype)
+    prd[0] = cs[0]*0
+    prd[1] = cs[0]
+    for i in range(1, len(cs)):
+        prd[i + 1] = cs[i]
+        prd[i - 1] += cs[i]*i
+    return prd
+
+
+def hermemul(c1, c2):
+    """
+    Multiply one Hermite series by another.
+
+    Returns the product of two Hermite series `c1` * `c2`.  The arguments
+    are sequences of coefficients, from lowest order "term" to highest,
+    e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Hermite series coefficients representing their product.
+
+    See Also
+    --------
+    hermeadd, hermesub, hermediv, hermepow
+
+    Notes
+    -----
+    In general, the (polynomial) product of two C-series results in terms
+    that are not in the Hermite polynomial basis set.  Thus, to express
+    the product as a Hermite series, it is necessary to "re-project" the
+    product onto said basis set, which may produce "un-intuitive" (but
+    correct) results; see Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermemul
+    >>> hermemul([1, 2, 3], [0, 1, 2])
+    array([ 14.,  15.,  28.,   7.,   6.])
+
+    """
+    # s1, s2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+
+    if len(c1) > len(c2):
+        cs = c2
+        xs = c1
+    else:
+        cs = c1
+        xs = c2
+
+    if len(cs) == 1:
+        c0 = cs[0]*xs
+        c1 = 0
+    elif len(cs) == 2:
+        c0 = cs[0]*xs
+        c1 = cs[1]*xs
+    else :
+        nd = len(cs)
+        c0 = cs[-2]*xs
+        c1 = cs[-1]*xs
+        for i in range(3, len(cs) + 1) :
+            tmp = c0
+            nd =  nd - 1
+            c0 = hermesub(cs[-i]*xs, c1*(nd - 1))
+            c1 = hermeadd(tmp, hermemulx(c1))
+    return hermeadd(c0, hermemulx(c1))
+
+
+def hermediv(c1, c2):
+    """
+    Divide one Hermite series by another.
+
+    Returns the quotient-with-remainder of two Hermite series
+    `c1` / `c2`.  The arguments are sequences of coefficients from lowest
+    order "term" to highest, e.g., [1,2,3] represents the series
+    ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Hermite series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    [quo, rem] : ndarrays
+        Of Hermite series coefficients representing the quotient and
+        remainder.
+
+    See Also
+    --------
+    hermeadd, hermesub, hermemul, hermepow
+
+    Notes
+    -----
+    In general, the (polynomial) division of one Hermite series by another
+    results in quotient and remainder terms that are not in the Hermite
+    polynomial basis set.  Thus, to express these results as a Hermite
+    series, it is necessary to "re-project" the results onto the Hermite
+    basis set, which may produce "un-intuitive" (but correct) results; see
+    Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermediv
+    >>> hermediv([ 14.,  15.,  28.,   7.,   6.], [0, 1, 2])
+    (array([ 1.,  2.,  3.]), array([ 0.]))
+    >>> hermediv([ 15.,  17.,  28.,   7.,   6.], [0, 1, 2])
+    (array([ 1.,  2.,  3.]), array([ 1.,  2.]))
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if c2[-1] == 0 :
+        raise ZeroDivisionError()
+
+    lc1 = len(c1)
+    lc2 = len(c2)
+    if lc1 < lc2 :
+        return c1[:1]*0, c1
+    elif lc2 == 1 :
+        return c1/c2[-1], c1[:1]*0
+    else :
+        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
+        rem = c1
+        for i in range(lc1 - lc2, - 1, -1):
+            p = hermemul([0]*i + [1], c2)
+            q = rem[-1]/p[-1]
+            rem = rem[:-1] - q*p[:-1]
+            quo[i] = q
+        return quo, pu.trimseq(rem)
+
+
+def hermepow(cs, pow, maxpower=16) :
+    """Raise a Hermite series to a power.
+
+    Returns the Hermite series `cs` raised to the power `pow`. The
+    arguement `cs` is a sequence of coefficients ordered from low to high.
+    i.e., [1,2,3] is the series  ``P_0 + 2*P_1 + 3*P_2.``
+
+    Parameters
+    ----------
+    cs : array_like
+        1d array of Hermite series coefficients ordered from low to
+        high.
+    pow : integer
+        Power to which the series will be raised
+    maxpower : integer, optional
+        Maximum power allowed. This is mainly to limit growth of the series
+        to umanageable size. Default is 16
+
+    Returns
+    -------
+    coef : ndarray
+        Hermite series of power.
+
+    See Also
+    --------
+    hermeadd, hermesub, hermemul, hermediv
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermepow
+    >>> hermepow([1, 2, 3], 2)
+    array([ 23.,  28.,  46.,  12.,   9.])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    power = int(pow)
+    if power != pow or power < 0 :
+        raise ValueError("Power must be a non-negative integer.")
+    elif maxpower is not None and power > maxpower :
+        raise ValueError("Power is too large")
+    elif power == 0 :
+        return np.array([1], dtype=cs.dtype)
+    elif power == 1 :
+        return cs
+    else :
+        # This can be made more efficient by using powers of two
+        # in the usual way.
+        prd = cs
+        for i in range(2, power + 1) :
+            prd = hermemul(prd, cs)
+        return prd
+
+
+def hermeder(cs, m=1, scl=1) :
+    """
+    Differentiate a Hermite series.
+
+    Returns the series `cs` differentiated `m` times.  At each iteration the
+    result is multiplied by `scl` (the scaling factor is for use in a linear
+    change of variable).  The argument `cs` is the sequence of coefficients
+    from lowest order "term" to highest, e.g., [1,2,3] represents the series
+    ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    cs: array_like
+        1-d array of Hermite series coefficients ordered from low to high.
+    m : int, optional
+        Number of derivatives taken, must be non-negative. (Default: 1)
+    scl : scalar, optional
+        Each differentiation is multiplied by `scl`.  The end result is
+        multiplication by ``scl**m``.  This is for use in a linear change of
+        variable. (Default: 1)
+
+    Returns
+    -------
+    der : ndarray
+        Hermite series of the derivative.
+
+    See Also
+    --------
+    hermeint
+
+    Notes
+    -----
+    In general, the result of differentiating a Hermite series does not
+    resemble the same operation on a power series. Thus the result of this
+    function may be "un-intuitive," albeit correct; see Examples section
+    below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermeder
+    >>> hermeder([ 1.,  1.,  1.,  1.])
+    array([ 1.,  2.,  3.])
+    >>> hermeder([-0.25,  1.,  1./2.,  1./3.,  1./4 ], m=2)
+    array([ 1.,  2.,  3.])
+
+    """
+    cnt = int(m)
+
+    if cnt != m:
+        raise ValueError, "The order of derivation must be integer"
+    if cnt < 0 :
+        raise ValueError, "The order of derivation must be non-negative"
+
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if cnt == 0:
+        return cs
+    elif cnt >= len(cs):
+        return cs[:1]*0
+    else :
+        for i in range(cnt):
+            n = len(cs) - 1
+            cs *= scl
+            der = np.empty(n, dtype=cs.dtype)
+            for j in range(n, 0, -1):
+                der[j - 1] = j*cs[j]
+            cs = der
+        return cs
+
+
+def hermeint(cs, m=1, k=[], lbnd=0, scl=1):
+    """
+    Integrate a Hermite series.
+
+    Returns a Hermite series that is the Hermite series `cs`, integrated
+    `m` times from `lbnd` to `x`.  At each iteration the resulting series
+    is **multiplied** by `scl` and an integration constant, `k`, is added.
+    The scaling factor is for use in a linear change of variable.  ("Buyer
+    beware": note that, depending on what one is doing, one may want `scl`
+    to be the reciprocal of what one might expect; for more information,
+    see the Notes section below.)  The argument `cs` is a sequence of
+    coefficients, from lowest order Hermite series "term" to highest,
+    e.g., [1,2,3] represents the series :math:`P_0(x) + 2P_1(x) + 3P_2(x)`.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Hermite series coefficients, ordered from low to high.
+    m : int, optional
+        Order of integration, must be positive. (Default: 1)
+    k : {[], list, scalar}, optional
+        Integration constant(s).  The value of the first integral at
+        ``lbnd`` is the first value in the list, the value of the second
+        integral at ``lbnd`` is the second value, etc.  If ``k == []`` (the
+        default), all constants are set to zero.  If ``m == 1``, a single
+        scalar can be given instead of a list.
+    lbnd : scalar, optional
+        The lower bound of the integral. (Default: 0)
+    scl : scalar, optional
+        Following each integration the result is *multiplied* by `scl`
+        before the integration constant is added. (Default: 1)
+
+    Returns
+    -------
+    S : ndarray
+        Hermite series coefficients of the integral.
+
+    Raises
+    ------
+    ValueError
+        If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or
+        ``np.isscalar(scl) == False``.
+
+    See Also
+    --------
+    hermeder
+
+    Notes
+    -----
+    Note that the result of each integration is *multiplied* by `scl`.
+    Why is this important to note?  Say one is making a linear change of
+    variable :math:`u = ax + b` in an integral relative to `x`.  Then
+    :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a`
+    - perhaps not what one would have first thought.
+
+    Also note that, in general, the result of integrating a C-series needs
+    to be "re-projected" onto the C-series basis set.  Thus, typically,
+    the result of this function is "un-intuitive," albeit correct; see
+    Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermeint
+    >>> hermeint([1, 2, 3]) # integrate once, value 0 at 0.
+    array([ 1.,  1.,  1.,  1.])
+    >>> hermeint([1, 2, 3], m=2) # integrate twice, value & deriv 0 at 0 
+    array([-0.25      ,  1.        ,  0.5       ,  0.33333333,  0.25      ])
+    >>> hermeint([1, 2, 3], k=1) # integrate once, value 1 at 0.
+    array([ 2.,  1.,  1.,  1.])
+    >>> hermeint([1, 2, 3], lbnd=-1) # integrate once, value 0 at -1
+    array([-1.,  1.,  1.,  1.])
+    >>> hermeint([1, 2, 3], m=2, k=[1,2], lbnd=-1)
+    array([ 1.83333333,  0.        ,  0.5       ,  0.33333333,  0.25      ])
+
+    """
+    cnt = int(m)
+    if np.isscalar(k) :
+        k = [k]
+
+    if cnt != m:
+        raise ValueError, "The order of integration must be integer"
+    if cnt < 0 :
+        raise ValueError, "The order of integration must be non-negative"
+    if len(k) > cnt :
+        raise ValueError, "Too many integration constants"
+
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if cnt == 0:
+        return cs
+
+    k = list(k) + [0]*(cnt - len(k))
+    for i in range(cnt) :
+        n = len(cs)
+        cs *= scl
+        if n == 1 and cs[0] == 0:
+            cs[0] += k[i]
+        else:
+            tmp = np.empty(n + 1, dtype=cs.dtype)
+            tmp[0] = cs[0]*0
+            tmp[1] = cs[0]
+            for j in range(1, n):
+                tmp[j + 1] = cs[j]/(j + 1)
+            tmp[0] += k[i] - hermeval(lbnd, tmp)
+            cs = tmp
+    return cs
+
+
+def hermeval(x, cs):
+    """Evaluate a Hermite series.
+
+    If `cs` is of length `n`, this function returns :
+
+    ``p(x) = cs[0]*P_0(x) + cs[1]*P_1(x) + ... + cs[n-1]*P_{n-1}(x)``
+
+    If x is a sequence or array then p(x) will have the same shape as x.
+    If r is a ring_like object that supports multiplication and addition
+    by the values in `cs`, then an object of the same type is returned.
+
+    Parameters
+    ----------
+    x : array_like, ring_like
+        Array of numbers or objects that support multiplication and
+        addition with themselves and with the elements of `cs`.
+    cs : array_like
+        1-d array of Hermite coefficients ordered from low to high.
+
+    Returns
+    -------
+    values : ndarray, ring_like
+        If the return is an ndarray then it has the same shape as `x`.
+
+    See Also
+    --------
+    hermefit
+
+    Examples
+    --------
+
+    Notes
+    -----
+    The evaluation uses Clenshaw recursion, aka synthetic division.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermeval
+    >>> coef = [1,2,3]
+    >>> hermeval(1, coef)
+    3.0
+    >>> hermeval([[1,2],[3,4]], coef)
+    array([[  3.,  14.],
+           [ 31.,  54.]])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if isinstance(x, tuple) or isinstance(x, list) :
+        x = np.asarray(x)
+
+    if len(cs) == 1 :
+        c0 = cs[0]
+        c1 = 0
+    elif len(cs) == 2 :
+        c0 = cs[0]
+        c1 = cs[1]
+    else :
+        nd = len(cs)
+        c0 = cs[-2]
+        c1 = cs[-1]
+        for i in range(3, len(cs) + 1) :
+            tmp = c0
+            nd =  nd - 1
+            c0 = cs[-i] - c1*(nd - 1)
+            c1 = tmp + c1*x
+    return c0 + c1*x
+
+
+def hermevander(x, deg) :
+    """Vandermonde matrix of given degree.
+
+    Returns the Vandermonde matrix of degree `deg` and sample points `x`.
+    This isn't a true Vandermonde matrix because `x` can be an arbitrary
+    ndarray and the Hermite polynomials aren't powers. If ``V`` is the
+    returned matrix and `x` is a 2d array, then the elements of ``V`` are
+    ``V[i,j,k] = P_k(x[i,j])``, where ``P_k`` is the Hermite polynomial
+    of degree ``k``.
+
+    Parameters
+    ----------
+    x : array_like
+        Array of points. The values are converted to double or complex
+        doubles. If x is scalar it is converted to a 1D array.
+    deg : integer
+        Degree of the resulting matrix.
+
+    Returns
+    -------
+    vander : Vandermonde matrix.
+        The shape of the returned matrix is ``x.shape + (deg+1,)``. The last
+        index is the degree.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermevander
+    >>> x = np.array([-1, 0, 1])
+    >>> hermevander(x, 3)
+    array([[ 1., -1.,  0.,  2.],
+           [ 1.,  0., -1., -0.],
+           [ 1.,  1.,  0., -2.]])
+
+    """
+    ideg = int(deg)
+    if ideg != deg:
+        raise ValueError("deg must be integer")
+    if ideg < 0:
+        raise ValueError("deg must be non-negative")
+
+    x = np.array(x, copy=0, ndmin=1) + 0.0
+    v = np.empty((ideg + 1,) + x.shape, dtype=x.dtype)
+    v[0] = x*0 + 1
+    if ideg > 0 :
+        v[1] = x
+        for i in range(2, ideg + 1) :
+            v[i] = (v[i-1]*x - v[i-2]*(i - 1))
+    return np.rollaxis(v, 0, v.ndim)
+
+
+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.
+
+    Parameters
+    ----------
+    x : array_like, shape (M,)
+        x-coordinates of the M sample points ``(x[i], y[i])``.
+    y : array_like, shape (M,) or (M, K)
+        y-coordinates of the sample points. Several data sets of sample
+        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
+    rcond : float, optional
+        Relative condition number of the fit. Singular values smaller than
+        this relative to the largest singular value will be ignored. The
+        default value is len(x)*eps, where eps is the relative precision of
+        the float type, about 2e-16 in most cases.
+    full : bool, optional
+        Switch determining nature of return value. When it is False (the
+        default) just the coefficients are returned, when True diagnostic
+        information from the singular value decomposition is also returned.
+    w : array_like, shape (`M`,), optional
+        Weights. If not None, the contribution of each point
+        ``(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.
+
+    Returns
+    -------
+    coef : ndarray, shape (M,) or (M, K)
+        Hermite coefficients ordered from low to high. If `y` was 2-D,
+        the coefficients for the data in column k  of `y` are in column
+        `k`.
+
+    [residuals, rank, singular_values, rcond] : present when `full` = True
+        Residuals of the least-squares fit, the effective rank of the
+        scaled Vandermonde matrix and its singular values, and the
+        specified value of `rcond`. For more details, see `linalg.lstsq`.
+
+    Warns
+    -----
+    RankWarning
+        The rank of the coefficient matrix in the least-squares fit is
+        deficient. The warning is only raised if `full` = False.  The
+        warnings can be turned off by
+
+        >>> import warnings
+        >>> warnings.simplefilter('ignore', RankWarning)
+
+    See Also
+    --------
+    hermeval : Evaluates a Hermite series.
+    hermevander : Vandermonde matrix of Hermite series.
+    polyfit : least squares fit using polynomials.
+    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 solution are the coefficients ``c[i]`` of the Hermite series
+    ``P(x)`` that minimizes the squared error
+
+    ``E = \\sum_j |y_j - P(x_j)|^2``.
+
+    This problem is solved by setting up as the overdetermined matrix
+    equation
+
+    ``V(x)*c = 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
+    decomposition of ``V``.
+
+    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.
+
+    References
+    ----------
+    .. [1] Wikipedia, "Curve fitting",
+           http://en.wikipedia.org/wiki/Curve_fitting
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermefit, hermeval
+    >>> x = np.linspace(-10, 10)
+    >>> err = np.random.randn(len(x))/10
+    >>> y = hermeval(x, [1, 2, 3]) + err
+    >>> hermefit(x, y, 2)
+    array([ 1.01690445,  1.99951418,  2.99948696])
+
+    """
+    order = int(deg) + 1
+    x = np.asarray(x) + 0.0
+    y = np.asarray(y) + 0.0
+
+    # check arguments.
+    if deg < 0 :
+        raise ValueError, "expected deg >= 0"
+    if x.ndim != 1:
+        raise TypeError, "expected 1D vector for x"
+    if x.size == 0:
+        raise TypeError, "expected non-empty vector for x"
+    if y.ndim < 1 or y.ndim > 2 :
+        raise TypeError, "expected 1D or 2D array for y"
+    if len(x) != len(y):
+        raise TypeError, "expected x and y to have same length"
+
+    # set up the least squares matrices
+    lhs = hermevander(x, deg)
+    rhs = y
+    if w is not None:
+        w = np.asarray(w) + 0.0
+        if w.ndim != 1:
+            raise TypeError, "expected 1D vector for w"
+        if len(x) != len(w):
+            raise TypeError, "expected x and w to have same length"
+        # apply weights
+        if rhs.ndim == 2:
+            lhs *= w[:, np.newaxis]
+            rhs *= w[:, np.newaxis]
+        else:
+            lhs *= w[:, np.newaxis]
+            rhs *= w
+
+    # set rcond
+    if rcond is None :
+        rcond = len(x)*np.finfo(x.dtype).eps
+
+    # scale the design matrix and solve the least squares equation
+    scl = np.sqrt((lhs*lhs).sum(0))
+    c, resids, rank, s = la.lstsq(lhs/scl, rhs, rcond)
+    c = (c.T/scl).T
+
+    # warn on rank reduction
+    if rank != order and not full:
+        msg = "The fit may be poorly conditioned"
+        warnings.warn(msg, pu.RankWarning)
+
+    if full :
+        return c, [resids, rank, s, rcond]
+    else :
+        return c
+
+
+def hermeroots(cs):
+    """
+    Compute the roots of a Hermite series.
+
+    Return the roots (a.k.a "zeros") of the Hermite 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.
+
+    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.
+
+    See Also
+    --------
+    polyroots
+    chebroots
+
+    Notes
+    -----
+    Algorithm(s) used:
+
+    Remember: because the Hermite series basis set is different from the
+    "standard" basis set, the results of this function *may* not be what
+    one is expecting.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.hermite_e import hermeroots, hermefromroots
+    >>> coef = hermefromroots([-1, 0, 1])
+    >>> coef
+    array([ 0.,  2.,  0.,  1.])
+    >>> hermeroots(coef)
+    array([-1.,  0.,  1.])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) <= 1 :
+        return np.array([], dtype=cs.dtype)
+    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] = 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
+
+
+#
+# HermiteE series class
+#
+
+exec polytemplate.substitute(name='HermiteE', nick='herme', domain='[-1,1]')
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
new file mode 100644
index 000000000..b6389bf63
--- /dev/null
+++ b/numpy/polynomial/laguerre.py
@@ -0,0 +1,1146 @@
+"""
+Objects for dealing with Laguerre series.
+
+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`).
+
+Constants
+---------
+- `lagdomain` -- Laguerre series default domain, [-1,1].
+- `lagzero` -- Laguerre series that evaluates identically to 0.
+- `lagone` -- Laguerre series that evaluates identically to 1.
+- `lagx` -- Laguerre series for the identity map, ``f(x) = x``.
+
+Arithmetic
+----------
+- `lagmulx` -- multiply a Laguerre series in ``P_i(x)`` by ``x``.
+- `lagadd` -- add two Laguerre series.
+- `lagsub` -- subtract one Laguerre series from another.
+- `lagmul` -- multiply two Laguerre series.
+- `lagdiv` -- divide one Laguerre series by another.
+- `lagval` -- evaluate a Laguerre series at given points.
+
+Calculus
+--------
+- `lagder` -- differentiate a Laguerre series.
+- `lagint` -- integrate a Laguerre series.
+
+Misc Functions
+--------------
+- `lagfromroots` -- create a Laguerre series with specified roots.
+- `lagroots` -- find the roots of a Laguerre series.
+- `lagvander` -- Vandermonde-like matrix for Laguerre polynomials.
+- `lagfit` -- least-squares fit returning a Laguerre series.
+- `lagtrim` -- trim leading coefficients from a Laguerre series.
+- `lagline` -- Laguerre series of given straight line.
+- `lag2poly` -- convert a Laguerre series to a polynomial.
+- `poly2lag` -- convert a polynomial to a Laguerre series.
+
+Classes
+-------
+- `Laguerre` -- A Laguerre series class.
+
+See also
+--------
+`numpy.polynomial`
+
+"""
+from __future__ import division
+
+__all__ = ['lagzero', 'lagone', 'lagx', 'lagdomain', 'lagline',
+        'lagadd', 'lagsub', 'lagmulx', 'lagmul', 'lagdiv', 'lagval',
+        'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots',
+        'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre']
+
+import numpy as np
+import numpy.linalg as la
+import polyutils as pu
+import warnings
+from polytemplate import polytemplate
+
+lagtrim = pu.trimcoef
+
+def poly2lag(pol) :
+    """
+    poly2lag(pol)
+
+    Convert a polynomial to a Laguerre series.
+
+    Convert an array representing the coefficients of a polynomial (relative
+    to the "standard" basis) ordered from lowest degree to highest, to an
+    array of the coefficients of the equivalent Laguerre series, ordered
+    from lowest to highest degree.
+
+    Parameters
+    ----------
+    pol : array_like
+        1-d array containing the polynomial coefficients
+
+    Returns
+    -------
+    cs : ndarray
+        1-d array containing the coefficients of the equivalent Laguerre
+        series.
+
+    See Also
+    --------
+    lag2poly
+
+    Notes
+    -----
+    The easy way to do conversions between polynomial basis sets
+    is to use the convert method of a class instance.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import poly2lag
+    >>> poly2lag(np.arange(4))
+    array([ 23., -63.,  58., -18.])
+
+    """
+    [pol] = pu.as_series([pol])
+    deg = len(pol) - 1
+    res = 0
+    for i in range(deg, -1, -1) :
+        res = lagadd(lagmulx(res), pol[i])
+    return res
+
+
+def lag2poly(cs) :
+    """
+    Convert a Laguerre series to a polynomial.
+
+    Convert an array representing the coefficients of a Laguerre series,
+    ordered from lowest degree to highest, to an array of the coefficients
+    of the equivalent polynomial (relative to the "standard" basis) ordered
+    from lowest to highest degree.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array containing the Laguerre series coefficients, ordered
+        from lowest order term to highest.
+
+    Returns
+    -------
+    pol : ndarray
+        1-d array containing the coefficients of the equivalent polynomial
+        (relative to the "standard" basis) ordered from lowest order term
+        to highest.
+
+    See Also
+    --------
+    poly2lag
+
+    Notes
+    -----
+    The easy way to do conversions between polynomial basis sets
+    is to use the convert method of a class instance.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lag2poly
+    >>> lag2poly([ 23., -63.,  58., -18.])
+    array([ 0.,  1.,  2.,  3.])
+
+    """
+    from polynomial import polyadd, polysub, polymulx
+
+    [cs] = pu.as_series([cs])
+    n = len(cs)
+    if n == 1:
+        return cs
+    else:
+        c0 = cs[-2]
+        c1 = cs[-1]
+        # i is the current degree of c1
+        for i in range(n - 1, 1, -1):
+            tmp = c0
+            c0 = polysub(cs[i - 2], (c1*(i - 1))/i)
+            c1 = polyadd(tmp, polysub((2*i - 1)*c1, polymulx(c1))/i)
+        return polyadd(c0, polysub(c1, polymulx(c1)))
+
+#
+# These are constant arrays are of integer type so as to be compatible
+# with the widest range of other types, such as Decimal.
+#
+
+# Laguerre
+lagdomain = np.array([0,1])
+
+# Laguerre coefficients representing zero.
+lagzero = np.array([0])
+
+# Laguerre coefficients representing one.
+lagone = np.array([1])
+
+# Laguerre coefficients representing the identity x.
+lagx = np.array([1, -1])
+
+
+def lagline(off, scl) :
+    """
+    Laguerre series whose graph is a straight line.
+
+
+
+    Parameters
+    ----------
+    off, scl : scalars
+        The specified line is given by ``off + scl*x``.
+
+    Returns
+    -------
+    y : ndarray
+        This module's representation of the Laguerre series for
+        ``off + scl*x``.
+
+    See Also
+    --------
+    polyline, chebline
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagline, lagval
+    >>> lagval(0,lagline(3, 2))
+    3.0
+    >>> lagval(1,lagline(3, 2))
+    5.0
+
+    """
+    if scl != 0 :
+        return np.array([off + scl, -scl])
+    else :
+        return np.array([off])
+
+
+def lagfromroots(roots) :
+    """
+    Generate a Laguerre series with the given roots.
+
+    Return the array of coefficients for the P-series whose roots (a.k.a.
+    "zeros") are given by *roots*.  The returned array of coefficients is
+    ordered from lowest order "term" to highest, and zeros of multiplicity
+    greater than one must be included in *roots* a number of times equal
+    to their multiplicity (e.g., if `2` is a root of multiplicity three,
+    then [2,2,2] must be in *roots*).
+
+    Parameters
+    ----------
+    roots : array_like
+        Sequence containing the roots.
+
+    Returns
+    -------
+    out : ndarray
+        1-d array of the Laguerre series coefficients, ordered from low to
+        high.  If all roots are real, ``out.dtype`` is a float type;
+        otherwise, ``out.dtype`` is a complex type, even if all the
+        coefficients in the result are real (see Examples below).
+
+    See Also
+    --------
+    polyfromroots, chebfromroots
+
+    Notes
+    -----
+    What is returned are the :math:`c_i` such that:
+
+    .. math::
+
+        \\sum_{i=0}^{n} c_i*P_i(x) = \\prod_{i=0}^{n} (x - roots[i])
+
+    where ``n == len(roots)`` and :math:`P_i(x)` is the `i`-th Laguerre
+    (basis) polynomial over the domain `[-1,1]`.  Note that, unlike
+    `polyfromroots`, due to the nature of the Laguerre basis set, the
+    above identity *does not* imply :math:`c_n = 1` identically (see
+    Examples).
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagfromroots, lagval
+    >>> coef = lagfromroots((-1, 0, 1))
+    >>> lagval((-1, 0, 1), coef)
+    array([ 0.,  0.,  0.])
+    >>> coef = lagfromroots((-1j, 1j))
+    >>> lagval((-1j, 1j), coef)
+    array([ 0.+0.j,  0.+0.j])
+
+    """
+    if len(roots) == 0 :
+        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
+
+
+def lagadd(c1, c2):
+    """
+    Add one Laguerre series to another.
+
+    Returns the sum of two Laguerre series `c1` + `c2`.  The arguments
+    are sequences of coefficients ordered from lowest order term to
+    highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Laguerre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the Laguerre series of their sum.
+
+    See Also
+    --------
+    lagsub, lagmul, lagdiv, lagpow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the sum of two Laguerre series
+    is a Laguerre series (without having to "reproject" the result onto
+    the basis set) so addition, just like that of "standard" polynomials,
+    is simply "component-wise."
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagadd
+    >>> lagadd([1, 2, 3], [1, 2, 3, 4])
+    array([ 2.,  4.,  6.,  4.])
+
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if len(c1) > len(c2) :
+        c1[:c2.size] += c2
+        ret = c1
+    else :
+        c2[:c1.size] += c1
+        ret = c2
+    return pu.trimseq(ret)
+
+
+def lagsub(c1, c2):
+    """
+    Subtract one Laguerre series from another.
+
+    Returns the difference of two Laguerre series `c1` - `c2`.  The
+    sequences of coefficients are from lowest order term to highest, i.e.,
+    [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Laguerre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Laguerre series coefficients representing their difference.
+
+    See Also
+    --------
+    lagadd, lagmul, lagdiv, lagpow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the difference of two Laguerre
+    series is a Laguerre series (without having to "reproject" the result
+    onto the basis set) so subtraction, just like that of "standard"
+    polynomials, is simply "component-wise."
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagsub
+    >>> lagsub([1, 2, 3, 4], [1, 2, 3])
+    array([ 0.,  0.,  0.,  4.])
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if len(c1) > len(c2) :
+        c1[:c2.size] -= c2
+        ret = c1
+    else :
+        c2 = -c2
+        c2[:c1.size] += c1
+        ret = c2
+    return pu.trimseq(ret)
+
+
+def lagmulx(cs):
+    """Multiply a Laguerre series by x.
+
+    Multiply the Laguerre series `cs` by x, where x is the independent
+    variable.
+
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Laguerre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the result of the multiplication.
+
+    Notes
+    -----
+    The multiplication uses the recursion relationship for Laguerre
+    polynomials in the form
+
+    .. math::
+
+    xP_i(x) = (-(i + 1)*P_{i + 1}(x) + (2i + 1)P_{i}(x) - iP_{i - 1}(x))
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagmulx
+    >>> lagmulx([1, 2, 3])
+    array([ -1.,  -1.,  11.,  -9.])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    # The zero series needs special treatment
+    if len(cs) == 1 and cs[0] == 0:
+        return cs
+
+    prd = np.empty(len(cs) + 1, dtype=cs.dtype)
+    prd[0] = cs[0]
+    prd[1] = -cs[0]
+    for i in range(1, len(cs)):
+        prd[i + 1] = -cs[i]*(i + 1)
+        prd[i] += cs[i]*(2*i + 1)
+        prd[i - 1] -= cs[i]*i
+    return prd
+
+
+def lagmul(c1, c2):
+    """
+    Multiply one Laguerre series by another.
+
+    Returns the product of two Laguerre series `c1` * `c2`.  The arguments
+    are sequences of coefficients, from lowest order "term" to highest,
+    e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Laguerre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Laguerre series coefficients representing their product.
+
+    See Also
+    --------
+    lagadd, lagsub, lagdiv, lagpow
+
+    Notes
+    -----
+    In general, the (polynomial) product of two C-series results in terms
+    that are not in the Laguerre polynomial basis set.  Thus, to express
+    the product as a Laguerre series, it is necessary to "re-project" the
+    product onto said basis set, which may produce "un-intuitive" (but
+    correct) results; see Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagmul
+    >>> lagmul([1, 2, 3], [0, 1, 2])
+    array([  8., -13.,  38., -51.,  36.])
+
+    """
+    # s1, s2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+
+    if len(c1) > len(c2):
+        cs = c2
+        xs = c1
+    else:
+        cs = c1
+        xs = c2
+
+    if len(cs) == 1:
+        c0 = cs[0]*xs
+        c1 = 0
+    elif len(cs) == 2:
+        c0 = cs[0]*xs
+        c1 = cs[1]*xs
+    else :
+        nd = len(cs)
+        c0 = cs[-2]*xs
+        c1 = cs[-1]*xs
+        for i in range(3, len(cs) + 1) :
+            tmp = c0
+            nd =  nd - 1
+            c0 = lagsub(cs[-i]*xs, (c1*(nd - 1))/nd)
+            c1 = lagadd(tmp, lagsub((2*nd - 1)*c1, lagmulx(c1))/nd)
+    return lagadd(c0, lagsub(c1, lagmulx(c1)))
+
+
+def lagdiv(c1, c2):
+    """
+    Divide one Laguerre series by another.
+
+    Returns the quotient-with-remainder of two Laguerre series
+    `c1` / `c2`.  The arguments are sequences of coefficients from lowest
+    order "term" to highest, e.g., [1,2,3] represents the series
+    ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    c1, c2 : array_like
+        1-d arrays of Laguerre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    [quo, rem] : ndarrays
+        Of Laguerre series coefficients representing the quotient and
+        remainder.
+
+    See Also
+    --------
+    lagadd, lagsub, lagmul, lagpow
+
+    Notes
+    -----
+    In general, the (polynomial) division of one Laguerre series by another
+    results in quotient and remainder terms that are not in the Laguerre
+    polynomial basis set.  Thus, to express these results as a Laguerre
+    series, it is necessary to "re-project" the results onto the Laguerre
+    basis set, which may produce "un-intuitive" (but correct) results; see
+    Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagdiv
+    >>> lagdiv([  8., -13.,  38., -51.,  36.], [0, 1, 2])
+    (array([ 1.,  2.,  3.]), array([ 0.]))
+    >>> lagdiv([  9., -12.,  38., -51.,  36.], [0, 1, 2])
+    (array([ 1.,  2.,  3.]), array([ 1.,  1.]))
+
+    """
+    # c1, c2 are trimmed copies
+    [c1, c2] = pu.as_series([c1, c2])
+    if c2[-1] == 0 :
+        raise ZeroDivisionError()
+
+    lc1 = len(c1)
+    lc2 = len(c2)
+    if lc1 < lc2 :
+        return c1[:1]*0, c1
+    elif lc2 == 1 :
+        return c1/c2[-1], c1[:1]*0
+    else :
+        quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
+        rem = c1
+        for i in range(lc1 - lc2, - 1, -1):
+            p = lagmul([0]*i + [1], c2)
+            q = rem[-1]/p[-1]
+            rem = rem[:-1] - q*p[:-1]
+            quo[i] = q
+        return quo, pu.trimseq(rem)
+
+
+def lagpow(cs, pow, maxpower=16) :
+    """Raise a Laguerre series to a power.
+
+    Returns the Laguerre series `cs` raised to the power `pow`. The
+    arguement `cs` is a sequence of coefficients ordered from low to high.
+    i.e., [1,2,3] is the series  ``P_0 + 2*P_1 + 3*P_2.``
+
+    Parameters
+    ----------
+    cs : array_like
+        1d array of Laguerre series coefficients ordered from low to
+        high.
+    pow : integer
+        Power to which the series will be raised
+    maxpower : integer, optional
+        Maximum power allowed. This is mainly to limit growth of the series
+        to umanageable size. Default is 16
+
+    Returns
+    -------
+    coef : ndarray
+        Laguerre series of power.
+
+    See Also
+    --------
+    lagadd, lagsub, lagmul, lagdiv
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagpow
+    >>> lagpow([1, 2, 3], 2)
+    array([ 14., -16.,  56., -72.,  54.])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    power = int(pow)
+    if power != pow or power < 0 :
+        raise ValueError("Power must be a non-negative integer.")
+    elif maxpower is not None and power > maxpower :
+        raise ValueError("Power is too large")
+    elif power == 0 :
+        return np.array([1], dtype=cs.dtype)
+    elif power == 1 :
+        return cs
+    else :
+        # This can be made more efficient by using powers of two
+        # in the usual way.
+        prd = cs
+        for i in range(2, power + 1) :
+            prd = lagmul(prd, cs)
+        return prd
+
+
+def lagder(cs, m=1, scl=1) :
+    """
+    Differentiate a Laguerre series.
+
+    Returns the series `cs` differentiated `m` times.  At each iteration the
+    result is multiplied by `scl` (the scaling factor is for use in a linear
+    change of variable).  The argument `cs` is the sequence of coefficients
+    from lowest order "term" to highest, e.g., [1,2,3] represents the series
+    ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    cs: array_like
+        1-d array of Laguerre series coefficients ordered from low to high.
+    m : int, optional
+        Number of derivatives taken, must be non-negative. (Default: 1)
+    scl : scalar, optional
+        Each differentiation is multiplied by `scl`.  The end result is
+        multiplication by ``scl**m``.  This is for use in a linear change of
+        variable. (Default: 1)
+
+    Returns
+    -------
+    der : ndarray
+        Laguerre series of the derivative.
+
+    See Also
+    --------
+    lagint
+
+    Notes
+    -----
+    In general, the result of differentiating a Laguerre series does not
+    resemble the same operation on a power series. Thus the result of this
+    function may be "un-intuitive," albeit correct; see Examples section
+    below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagder
+    >>> lagder([ 1.,  1.,  1., -3.])
+    array([ 1.,  2.,  3.])
+    >>> lagder([ 1.,  0.,  0., -4.,  3.], m=2)
+    array([ 1.,  2.,  3.])
+
+    """
+    cnt = int(m)
+
+    if cnt != m:
+        raise ValueError, "The order of derivation must be integer"
+    if cnt < 0 :
+        raise ValueError, "The order of derivation must be non-negative"
+
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if cnt == 0:
+        return cs
+    elif cnt >= len(cs):
+        return cs[:1]*0
+    else :
+        for i in range(cnt):
+            n = len(cs) - 1
+            cs *= scl
+            der = np.empty(n, dtype=cs.dtype)
+            for j in range(n, 0, -1):
+                der[j - 1] = -cs[j]
+                cs[j - 1] += cs[j]
+            cs = der
+        return cs
+
+
+def lagint(cs, m=1, k=[], lbnd=0, scl=1):
+    """
+    Integrate a Laguerre series.
+
+    Returns a Laguerre series that is the Laguerre series `cs`, integrated
+    `m` times from `lbnd` to `x`.  At each iteration the resulting series
+    is **multiplied** by `scl` and an integration constant, `k`, is added.
+    The scaling factor is for use in a linear change of variable.  ("Buyer
+    beware": note that, depending on what one is doing, one may want `scl`
+    to be the reciprocal of what one might expect; for more information,
+    see the Notes section below.)  The argument `cs` is a sequence of
+    coefficients, from lowest order Laguerre series "term" to highest,
+    e.g., [1,2,3] represents the series :math:`P_0(x) + 2P_1(x) + 3P_2(x)`.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Laguerre series coefficients, ordered from low to high.
+    m : int, optional
+        Order of integration, must be positive. (Default: 1)
+    k : {[], list, scalar}, optional
+        Integration constant(s).  The value of the first integral at
+        ``lbnd`` is the first value in the list, the value of the second
+        integral at ``lbnd`` is the second value, etc.  If ``k == []`` (the
+        default), all constants are set to zero.  If ``m == 1``, a single
+        scalar can be given instead of a list.
+    lbnd : scalar, optional
+        The lower bound of the integral. (Default: 0)
+    scl : scalar, optional
+        Following each integration the result is *multiplied* by `scl`
+        before the integration constant is added. (Default: 1)
+
+    Returns
+    -------
+    S : ndarray
+        Laguerre series coefficients of the integral.
+
+    Raises
+    ------
+    ValueError
+        If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or
+        ``np.isscalar(scl) == False``.
+
+    See Also
+    --------
+    lagder
+
+    Notes
+    -----
+    Note that the result of each integration is *multiplied* by `scl`.
+    Why is this important to note?  Say one is making a linear change of
+    variable :math:`u = ax + b` in an integral relative to `x`.  Then
+    :math:`dx = du/a`, so one will need to set `scl` equal to :math:`1/a`
+    - perhaps not what one would have first thought.
+
+    Also note that, in general, the result of integrating a C-series needs
+    to be "re-projected" onto the C-series basis set.  Thus, typically,
+    the result of this function is "un-intuitive," albeit correct; see
+    Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagint
+    >>> lagint([1,2,3])
+    array([ 1.,  1.,  1., -3.])
+    >>> lagint([1,2,3], m=2)
+    array([ 1.,  0.,  0., -4.,  3.])
+    >>> lagint([1,2,3], k=1)
+    array([ 2.,  1.,  1., -3.])
+    >>> lagint([1,2,3], lbnd=-1)
+    array([ 11.5,   1. ,   1. ,  -3. ])
+    >>> lagint([1,2], m=2, k=[1,2], lbnd=-1)
+    array([ 11.16666667,  -5.        ,  -3.        ,   2.        ])
+
+    """
+    cnt = int(m)
+    if np.isscalar(k) :
+        k = [k]
+
+    if cnt != m:
+        raise ValueError, "The order of integration must be integer"
+    if cnt < 0 :
+        raise ValueError, "The order of integration must be non-negative"
+    if len(k) > cnt :
+        raise ValueError, "Too many integration constants"
+
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if cnt == 0:
+        return cs
+
+    k = list(k) + [0]*(cnt - len(k))
+    for i in range(cnt) :
+        n = len(cs)
+        cs *= scl
+        if n == 1 and cs[0] == 0:
+            cs[0] += k[i]
+        else:
+            tmp = np.empty(n + 1, dtype=cs.dtype)
+            tmp[0] = cs[0]
+            tmp[1] = -cs[0]
+            for j in range(1, n):
+                tmp[j] += cs[j]
+                tmp[j + 1] = -cs[j]
+            tmp[0] += k[i] - lagval(lbnd, tmp)
+            cs = tmp
+    return cs
+
+
+def lagval(x, cs):
+    """Evaluate a Laguerre series.
+
+    If `cs` is of length `n`, this function returns :
+
+    ``p(x) = cs[0]*P_0(x) + cs[1]*P_1(x) + ... + cs[n-1]*P_{n-1}(x)``
+
+    If x is a sequence or array then p(x) will have the same shape as x.
+    If r is a ring_like object that supports multiplication and addition
+    by the values in `cs`, then an object of the same type is returned.
+
+    Parameters
+    ----------
+    x : array_like, ring_like
+        Array of numbers or objects that support multiplication and
+        addition with themselves and with the elements of `cs`.
+    cs : array_like
+        1-d array of Laguerre coefficients ordered from low to high.
+
+    Returns
+    -------
+    values : ndarray, ring_like
+        If the return is an ndarray then it has the same shape as `x`.
+
+    See Also
+    --------
+    lagfit
+
+    Examples
+    --------
+
+    Notes
+    -----
+    The evaluation uses Clenshaw recursion, aka synthetic division.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagval
+    >>> coef = [1,2,3]
+    >>> lagval(1, coef)
+    -0.5
+    >>> lagval([[1,2],[3,4]], coef)
+    array([[-0.5, -4. ],
+           [-4.5, -2. ]])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if isinstance(x, tuple) or isinstance(x, list) :
+        x = np.asarray(x)
+
+    if len(cs) == 1 :
+        c0 = cs[0]
+        c1 = 0
+    elif len(cs) == 2 :
+        c0 = cs[0]
+        c1 = cs[1]
+    else :
+        nd = len(cs)
+        c0 = cs[-2]
+        c1 = cs[-1]
+        for i in range(3, len(cs) + 1) :
+            tmp = c0
+            nd =  nd - 1
+            c0 = cs[-i] - (c1*(nd - 1))/nd
+            c1 = tmp + (c1*((2*nd - 1) - x))/nd
+    return c0 + c1*(1 - x)
+
+
+def lagvander(x, deg) :
+    """Vandermonde matrix of given degree.
+
+    Returns the Vandermonde matrix of degree `deg` and sample points `x`.
+    This isn't a true Vandermonde matrix because `x` can be an arbitrary
+    ndarray and the Laguerre polynomials aren't powers. If ``V`` is the
+    returned matrix and `x` is a 2d array, then the elements of ``V`` are
+    ``V[i,j,k] = P_k(x[i,j])``, where ``P_k`` is the Laguerre polynomial
+    of degree ``k``.
+
+    Parameters
+    ----------
+    x : array_like
+        Array of points. The values are converted to double or complex
+        doubles. If x is scalar it is converted to a 1D array.
+    deg : integer
+        Degree of the resulting matrix.
+
+    Returns
+    -------
+    vander : Vandermonde matrix.
+        The shape of the returned matrix is ``x.shape + (deg+1,)``. The last
+        index is the degree.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagvander
+    >>> x = np.array([0, 1, 2])
+    >>> lagvander(x, 3)
+    array([[ 1.        ,  1.        ,  1.        ,  1.        ],
+           [ 1.        ,  0.        , -0.5       , -0.66666667],
+           [ 1.        , -1.        , -1.        , -0.33333333]])
+
+    """
+    ideg = int(deg)
+    if ideg != deg:
+        raise ValueError("deg must be integer")
+    if ideg < 0:
+        raise ValueError("deg must be non-negative")
+
+    x = np.array(x, copy=0, ndmin=1) + 0.0
+    v = np.empty((ideg + 1,) + x.shape, dtype=x.dtype)
+    v[0] = x*0 + 1
+    if ideg > 0 :
+        v[1] = 1 - x
+        for i in range(2, ideg + 1) :
+            v[i] = (v[i-1]*(2*i - 1 - x) - v[i-2]*(i - 1))/i
+    return np.rollaxis(v, 0, v.ndim)
+
+
+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.
+
+    Parameters
+    ----------
+    x : array_like, shape (M,)
+        x-coordinates of the M sample points ``(x[i], y[i])``.
+    y : array_like, shape (M,) or (M, K)
+        y-coordinates of the sample points. Several data sets of sample
+        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
+    rcond : float, optional
+        Relative condition number of the fit. Singular values smaller than
+        this relative to the largest singular value will be ignored. The
+        default value is len(x)*eps, where eps is the relative precision of
+        the float type, about 2e-16 in most cases.
+    full : bool, optional
+        Switch determining nature of return value. When it is False (the
+        default) just the coefficients are returned, when True diagnostic
+        information from the singular value decomposition is also returned.
+    w : array_like, shape (`M`,), optional
+        Weights. If not None, the contribution of each point
+        ``(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.
+
+    Returns
+    -------
+    coef : ndarray, shape (M,) or (M, K)
+        Laguerre coefficients ordered from low to high. If `y` was 2-D,
+        the coefficients for the data in column k  of `y` are in column
+        `k`.
+
+    [residuals, rank, singular_values, rcond] : present when `full` = True
+        Residuals of the least-squares fit, the effective rank of the
+        scaled Vandermonde matrix and its singular values, and the
+        specified value of `rcond`. For more details, see `linalg.lstsq`.
+
+    Warns
+    -----
+    RankWarning
+        The rank of the coefficient matrix in the least-squares fit is
+        deficient. The warning is only raised if `full` = False.  The
+        warnings can be turned off by
+
+        >>> import warnings
+        >>> warnings.simplefilter('ignore', RankWarning)
+
+    See Also
+    --------
+    lagval : Evaluates a Laguerre series.
+    lagvander : Vandermonde matrix of Laguerre series.
+    polyfit : least squares fit using polynomials.
+    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 solution are the coefficients ``c[i]`` of the Laguerre series
+    ``P(x)`` that minimizes the squared error
+
+    ``E = \\sum_j |y_j - P(x_j)|^2``.
+
+    This problem is solved by setting up as the overdetermined matrix
+    equation
+
+    ``V(x)*c = 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
+    decomposition of ``V``.
+
+    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.
+
+    References
+    ----------
+    .. [1] Wikipedia, "Curve fitting",
+           http://en.wikipedia.org/wiki/Curve_fitting
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagfit, lagval
+    >>> x = np.linspace(0, 10)
+    >>> err = np.random.randn(len(x))/10
+    >>> y = lagval(x, [1, 2, 3]) + err
+    >>> lagfit(x, y, 2)
+    array([ 0.96971004,  2.00193749,  3.00288744])
+
+    """
+    order = int(deg) + 1
+    x = np.asarray(x) + 0.0
+    y = np.asarray(y) + 0.0
+
+    # check arguments.
+    if deg < 0 :
+        raise ValueError, "expected deg >= 0"
+    if x.ndim != 1:
+        raise TypeError, "expected 1D vector for x"
+    if x.size == 0:
+        raise TypeError, "expected non-empty vector for x"
+    if y.ndim < 1 or y.ndim > 2 :
+        raise TypeError, "expected 1D or 2D array for y"
+    if len(x) != len(y):
+        raise TypeError, "expected x and y to have same length"
+
+    # set up the least squares matrices
+    lhs = lagvander(x, deg)
+    rhs = y
+    if w is not None:
+        w = np.asarray(w) + 0.0
+        if w.ndim != 1:
+            raise TypeError, "expected 1D vector for w"
+        if len(x) != len(w):
+            raise TypeError, "expected x and w to have same length"
+        # apply weights
+        if rhs.ndim == 2:
+            lhs *= w[:, np.newaxis]
+            rhs *= w[:, np.newaxis]
+        else:
+            lhs *= w[:, np.newaxis]
+            rhs *= w
+
+    # set rcond
+    if rcond is None :
+        rcond = len(x)*np.finfo(x.dtype).eps
+
+    # scale the design matrix and solve the least squares equation
+    scl = np.sqrt((lhs*lhs).sum(0))
+    c, resids, rank, s = la.lstsq(lhs/scl, rhs, rcond)
+    c = (c.T/scl).T
+
+    # warn on rank reduction
+    if rank != order and not full:
+        msg = "The fit may be poorly conditioned"
+        warnings.warn(msg, pu.RankWarning)
+
+    if full :
+        return c, [resids, rank, s, rcond]
+    else :
+        return c
+
+
+def lagroots(cs):
+    """
+    Compute the roots of a Laguerre series.
+
+    Return the roots (a.k.a "zeros") of the Laguerre 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 Laguerre series coefficients ordered from low to high.
+
+    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.
+
+    See Also
+    --------
+    polyroots
+    chebroots
+
+    Notes
+    -----
+    Algorithm(s) used:
+
+    Remember: because the Laguerre series basis set is different from the
+    "standard" basis set, the results of this function *may* not be what
+    one is expecting.
+
+    Examples
+    --------
+    >>> from numpy.polynomial.laguerre import lagroots, lagfromroots
+    >>> coef = lagfromroots([0, 1, 2])
+    >>> coef
+    array([  2.,  -8.,  12.,  -6.])
+    >>> lagroots(coef)
+    array([ -4.44089210e-16,   1.00000000e+00,   2.00000000e+00])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) <= 1 :
+        return np.array([], dtype=cs.dtype)
+    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
+
+
+#
+# Laguerre series class
+#
+
+exec polytemplate.substitute(name='Laguerre', nick='lag', domain='[-1,1]')
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 9aec256cd..3ea68a5e6 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -227,52 +227,6 @@ def legline(off, scl) :
         return np.array([off])
 
 
-def legtimesx(cs):
-    """Multiply a Legendre series by x.
-
-    Multiply the Legendre series `cs` by x, where x is the independent
-    variable.
-
-
-    Parameters
-    ----------
-    cs : array_like
-        1-d array of Legendre series coefficients ordered from low to
-        high.
-
-    Returns
-    -------
-    out : ndarray
-        Array representing the result of the multiplication.
-
-    Notes
-    -----
-    The multiplication uses the recursion relationship for Legendre
-    polynomials in the form
-
-    .. math::
-
-    xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1)
-
-    """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    # The zero series needs special treatment
-    if len(cs) == 1 and cs[0] == 0:
-        return cs
-
-    prd = np.empty(len(cs) + 1, dtype=cs.dtype)
-    prd[0] = cs[0]*0
-    prd[1] = cs[0]
-    for i in range(1, len(cs)):
-        j = i + 1
-        k = i - 1
-        s = i + j
-        prd[j] = (cs[i]*j)/s
-        prd[k] += (cs[i]*i)/s
-    return prd
-
-
 def legfromroots(roots) :
     """
     Generate a Legendre series with the given roots.
diff --git a/numpy/polynomial/polytemplate.py b/numpy/polynomial/polytemplate.py
index 2106ad84e..657b48508 100644
--- a/numpy/polynomial/polytemplate.py
+++ b/numpy/polynomial/polytemplate.py
@@ -31,16 +31,22 @@ class $name(pu.PolyBase) :
         $name coefficients, in increasing order.  For example,
         ``(1, 2, 3)`` implies ``P_0 + 2P_1 + 3P_2`` where the
         ``P_i`` are a graded polynomial basis.
-    domain : (2,) array_like
+    domain : (2,) array_like, optional
         Domain to use. The interval ``[domain[0], domain[1]]`` is mapped to
-        the interval ``$domain`` by shifting and scaling.
+        the interval ``[window[0], window[1]]`` by shifting and scaling.
+        The default value is $domain.
+    window : (2,) array_like, optional
+        Window, see ``domain`` for its use. The default value is $domain.
+        .. versionadded:: 1.6.0
 
     Attributes
     ----------
     coef : (N,) array
         $name coefficients, from low to high.
-    domain : (2,) array_like
-        Domain that is mapped to ``$domain``.
+    domain : (2,) array
+        Domain that is mapped to ``window``.
+    window : (2,) array
+        Window that ``domain`` is mapped to.
 
     Class Attributes
     ----------------
@@ -49,6 +55,8 @@ class $name(pu.PolyBase) :
         ``p(x)**n`` is allowed. This is to limit runaway polynomial size.
     domain : (2,) ndarray
         Default domain of the class.
+    window : (2,) ndarray
+        Default window of the class.
 
     Notes
     -----
@@ -65,25 +73,115 @@ class $name(pu.PolyBase) :
     maxpower = 16
     # Default domain
     domain = np.array($domain)
+    # Default window
+    window = np.array($domain)
     # Don't let participate in array operations. Value doesn't matter.
     __array_priority__ = 0
 
-    def __init__(self, coef, domain=$domain) :
-        [coef, domain] = pu.as_series([coef, domain], trim=False)
-        if len(domain) != 2 :
+    def has_samecoef(self, other):
+        """Check if coefficients match.
+
+        Parameters
+        ----------
+        other : class instance
+            The other class must have the ``coef`` attribute.
+
+        Returns
+        -------
+        bool : boolean
+            True if the coefficients are the same, False otherwise.
+
+        Notes
+        -----
+        .. versionadded:: 1.6.0
+
+        """
+        if len(self.coef) != len(other.coef):
+            return False
+        elif not np.all(self.coef == other.coef):
+            return False
+        else:
+            return True
+
+    def has_samedomain(self, other):
+        """Check if domains match.
+
+        Parameters
+        ----------
+        other : class instance
+            The other class must have the ``domain`` attribute.
+
+        Returns
+        -------
+        bool : boolean
+            True if the domains are the same, False otherwise.
+
+        Notes
+        -----
+        .. versionadded:: 1.6.0
+
+        """
+        return np.all(self.domain == other.domain)
+
+    def has_samewindow(self, other):
+        """Check if windows match.
+
+        Parameters
+        ----------
+        other : class instance
+            The other class must have the ``window`` attribute.
+
+        Returns
+        -------
+        bool : boolean
+            True if the windows are the same, False otherwise.
+
+        Notes
+        -----
+        .. versionadded:: 1.6.0
+
+        """
+        return np.all(self.window == other.window)
+
+    def has_samewindow(self, other):
+        """Check if windows match.
+
+        Parameters
+        ----------
+        other : class instance
+            The other class must have the ``window`` attribute.
+
+        Returns
+        -------
+        bool : boolean
+            True if the windows are the same, False otherwise.
+
+        """
+        return np.all(self.window == other.window)
+
+    def __init__(self, coef, domain=$domain, window=$domain) :
+        [coef, dom, win] = pu.as_series([coef, domain, window], trim=False)
+        if len(dom) != 2 :
             raise ValueError("Domain has wrong number of elements.")
+        if len(win) != 2 :
+            raise ValueError("Window has wrong number of elements.")
         self.coef = coef
-        self.domain = domain
+        self.domain = dom
+        self.window = win
 
     def __repr__(self):
-        format = "%s(%s, %s)"
+        format = "%s(%s, %s, %s)"
         coef = repr(self.coef)[6:-1]
         domain = repr(self.domain)[6:-1]
-        return format % ('$name', coef, domain)
+        window = repr(self.domain)[6:-1]
+        return format % ('$name', coef, domain, window)
 
     def __str__(self) :
-        format = "%s(%s, %s)"
-        return format % ('$nick', str(self.coef), str(self.domain))
+        format = "%s(%s, %s, %s)"
+        coef = str(self.coef)[6:-1]
+        domain = str(self.domain)[6:-1]
+        window = str(self.domain)[6:-1]
+        return format % ('$nick', coef, domain, window)
 
     # Pickle and copy
 
@@ -91,6 +189,7 @@ class $name(pu.PolyBase) :
         ret = self.__dict__.copy()
         ret['coef'] = self.coef.copy()
         ret['domain'] = self.domain.copy()
+        ret['window'] = self.window.copy()
         return ret
 
     def __setstate__(self, dict) :
@@ -99,11 +198,10 @@ class $name(pu.PolyBase) :
     # Call
 
     def __call__(self, arg) :
-        off, scl = pu.mapparms(self.domain, $domain)
+        off, scl = pu.mapparms(self.domain, self.window)
         arg = off + scl*arg
         return ${nick}val(arg, self.coef)
 
-
     def __iter__(self) :
         return iter(self.coef)
 
@@ -112,9 +210,8 @@ class $name(pu.PolyBase) :
 
     # Numeric properties.
 
-
     def __neg__(self) :
-        return self.__class__(-self.coef, self.domain)
+        return self.__class__(-self.coef, self.domain, self.window)
 
     def __pos__(self) :
         return self
@@ -122,7 +219,7 @@ class $name(pu.PolyBase) :
     def __add__(self, other) :
         """Returns sum"""
         if isinstance(other, self.__class__) :
-            if np.all(self.domain == other.domain) :
+            if self.has_samedomain(other) and self.has_samewindow(other):
                 coef = ${nick}add(self.coef, other.coef)
             else :
                 raise PolyDomainError()
@@ -131,12 +228,12 @@ class $name(pu.PolyBase) :
                 coef = ${nick}add(self.coef, other)
             except :
                 return NotImplemented
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __sub__(self, other) :
         """Returns difference"""
         if isinstance(other, self.__class__) :
-            if np.all(self.domain == other.domain) :
+            if self.has_samedomain(other) and self.has_samewindow(other):
                 coef = ${nick}sub(self.coef, other.coef)
             else :
                 raise PolyDomainError()
@@ -145,12 +242,12 @@ class $name(pu.PolyBase) :
                 coef = ${nick}sub(self.coef, other)
             except :
                 return NotImplemented
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __mul__(self, other) :
         """Returns product"""
         if isinstance(other, self.__class__) :
-            if np.all(self.domain == other.domain) :
+            if self.has_samedomain(other) and self.has_samewindow(other):
                 coef = ${nick}mul(self.coef, other.coef)
             else :
                 raise PolyDomainError()
@@ -159,7 +256,7 @@ class $name(pu.PolyBase) :
                 coef = ${nick}mul(self.coef, other)
             except :
                 return NotImplemented
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __div__(self, other):
         # set to __floordiv__ /.
@@ -175,10 +272,11 @@ class $name(pu.PolyBase) :
             else :
                 return NotImplemented
         elif np.isscalar(other) :
+            # this might be overly restrictive
             coef = self.coef/other
         else :
             return NotImplemented
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __floordiv__(self, other) :
         """Returns the quotient."""
@@ -192,12 +290,12 @@ class $name(pu.PolyBase) :
                 quo, rem = ${nick}div(self.coef, other)
             except :
                 return NotImplemented
-        return self.__class__(quo, self.domain)
+        return self.__class__(quo, self.domain, self.window)
 
     def __mod__(self, other) :
         """Returns the remainder."""
         if isinstance(other, self.__class__) :
-            if np.all(self.domain == other.domain) :
+            if self.has_samedomain(other) and self.has_samewindow(other):
                 quo, rem = ${nick}div(self.coef, other.coef)
             else :
                 raise PolyDomainError()
@@ -206,12 +304,12 @@ class $name(pu.PolyBase) :
                 quo, rem = ${nick}div(self.coef, other)
             except :
                 return NotImplemented
-        return self.__class__(rem, self.domain)
+        return self.__class__(rem, self.domain, self.window)
 
     def __divmod__(self, other) :
         """Returns quo, remainder"""
         if isinstance(other, self.__class__) :
-            if np.all(self.domain == other.domain) :
+            if self.has_samedomain(other) and self.has_samewindow(other):
                 quo, rem = ${nick}div(self.coef, other.coef)
             else :
                 raise PolyDomainError()
@@ -220,35 +318,37 @@ class $name(pu.PolyBase) :
                 quo, rem = ${nick}div(self.coef, other)
             except :
                 return NotImplemented
-        return self.__class__(quo, self.domain), self.__class__(rem, self.domain)
+        quo = self.__class__(quo, self.domain, self.window)
+        rem = self.__class__(rem, self.domain, self.window)
+        return quo, rem
 
     def __pow__(self, other) :
         try :
             coef = ${nick}pow(self.coef, other, maxpower = self.maxpower)
         except :
             raise
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __radd__(self, other) :
         try :
             coef = ${nick}add(other, self.coef)
         except :
             return NotImplemented
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __rsub__(self, other):
         try :
             coef = ${nick}sub(other, self.coef)
         except :
             return NotImplemented
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __rmul__(self, other) :
         try :
             coef = ${nick}mul(other, self.coef)
         except :
             return NotImplemented
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def __rdiv__(self, other):
         # set to __floordiv__ /.
@@ -263,46 +363,62 @@ class $name(pu.PolyBase) :
                 quo, rem = ${nick}div(other, self.coef[0])
             except :
                 return NotImplemented
-        return self.__class__(quo, self.domain)
+        return self.__class__(quo, self.domain, self.window)
 
     def __rfloordiv__(self, other) :
         try :
             quo, rem = ${nick}div(other, self.coef)
         except :
             return NotImplemented
-        return self.__class__(quo, self.domain)
+        return self.__class__(quo, self.domain, self.window)
 
     def __rmod__(self, other) :
         try :
             quo, rem = ${nick}div(other, self.coef)
         except :
             return NotImplemented
-        return self.__class__(rem, self.domain)
+        return self.__class__(rem, self.domain, self.window)
 
     def __rdivmod__(self, other) :
         try :
             quo, rem = ${nick}div(other, self.coef)
         except :
             return NotImplemented
-        return self.__class__(quo, self.domain), self.__class__(rem, self.domain)
+        quo = self.__class__(quo, self.domain, self.window)
+        rem = self.__class__(rem, self.domain, self.window)
+        return quo, rem
 
     # Enhance me
     # some augmented arithmetic operations could be added here
 
     def __eq__(self, other) :
         res = isinstance(other, self.__class__) \
-                and len(self.coef) == len(other.coef) \
-                and np.all(self.domain == other.domain) \
-                and np.all(self.coef == other.coef)
+                and self.has_samecoef(other) \
+                and self.has_samedomain(other) \
+                and self.has_samewindow(other)
         return res
 
     def __ne__(self, other) :
         return not self.__eq__(other)
 
     #
-    # Extra numeric functions.
+    # Extra methods.
     #
 
+    def copy(self) :
+        """Return a copy.
+
+        A new instance of $name is returned that has the same
+        coefficients and domain as the current instance.
+
+        Returns
+        -------
+        new_instance : $name
+            New instance of $name with the same coefficients and domain.
+
+        """
+        return self.__class__(self.coef, self.domain, self.window)
+
     def degree(self) :
         """The degree of the series.
 
@@ -340,68 +456,6 @@ class $name(pu.PolyBase) :
         """
         return self.truncate(deg + 1)
 
-    def convert(self, domain=None, kind=None) :
-        """Convert to different class and/or domain.
-
-        Parameters
-        ----------
-        domain : array_like, optional
-            The domain of the new series type instance. If the value is None,
-            then the default domain of `kind` is used.
-        kind : class, optional
-            The polynomial series type class to which the current instance
-            should be converted. If kind is None, then the class of the
-            current instance is used.
-
-        Returns
-        -------
-        new_series_instance : `kind`
-            The returned class can be of different type than the current
-            instance and/or have a different domain.
-
-        Notes
-        -----
-        Conversion between domains and class types can result in
-        numerically ill defined series.
-
-        Examples
-        --------
-
-        """
-        if kind is None :
-            kind = $name
-        if domain is None :
-            domain = kind.domain
-        return self(kind.identity(domain))
-
-    def mapparms(self) :
-        """Return the mapping parameters.
-
-        The returned values define a linear map ``off + scl*x`` that is
-        applied to the input arguments before the series is evaluated. The
-        of the map depend on the domain; if the current domain is equal to
-        the default domain ``$domain`` the resulting map is the identity.
-        If the coeffients of the ``$name`` instance are to be used
-        separately, then the linear function must be substituted for the
-        ``x`` in the standard representation of the base polynomials.
-
-        Returns
-        -------
-        off, scl : floats or complex
-            The mapping function is defined by ``off + scl*x``.
-
-        Notes
-        -----
-        If the current domain is the interval ``[l_1, r_1]`` and the default
-        interval is ``[l_2, r_2]``, then the linear mapping function ``L`` is
-        defined by the equations::
-
-            L(l_1) = l_2
-            L(r_1) = r_2
-
-        """
-        return pu.mapparms(self.domain, $domain)
-
     def trim(self, tol=0) :
         """Remove small leading coefficients
 
@@ -422,7 +476,8 @@ class $name(pu.PolyBase) :
             Contains the new set of coefficients.
 
         """
-        return self.__class__(pu.trimcoef(self.coef, tol), self.domain)
+        coef = pu.trimcoef(self.coef, tol)
+        return self.__class__(coef, self.domain, self.window)
 
     def truncate(self, size) :
         """Truncate series to length `size`.
@@ -448,23 +503,75 @@ class $name(pu.PolyBase) :
         if isize != size or isize < 1 :
             raise ValueError("size must be a positive integer")
         if isize >= len(self.coef) :
-            return self.__class__(self.coef, self.domain)
+            coef = self.coef
         else :
-            return self.__class__(self.coef[:isize], self.domain)
+            coef = self.coef[:isize]
+        return self.__class__(coef, self.domain, self.window)
 
-    def copy(self) :
-        """Return a copy.
+    def convert(self, domain=None, kind=None, window=None) :
+        """Convert to different class and/or domain.
 
-        A new instance of $name is returned that has the same
-        coefficients and domain as the current instance.
+        Parameters
+        ----------
+        domain : array_like, optional
+            The domain of the new series type instance. If the value is None,
+            then the default domain of `kind` is used.
+        kind : class, optional
+            The polynomial series type class to which the current instance
+            should be converted. If kind is None, then the class of the
+            current instance is used.
 
         Returns
         -------
-        new_instance : $name
-            New instance of $name with the same coefficients and domain.
+        new_series_instance : `kind`
+            The returned class can be of different type than the current
+            instance and/or have a different domain.
+
+        Notes
+        -----
+        Conversion between domains and class types can result in
+        numerically ill defined series.
+
+        Examples
+        --------
+
+        """
+        if kind is None:
+            kind = $name
+        if domain is None:
+            domain = kind.domain
+        if window is None:
+            window = kind.window
+        return self(kind.identity(domain, window=window))
+
+    def mapparms(self) :
+        """Return the mapping parameters.
+
+        The returned values define a linear map ``off + scl*x`` that is
+        applied to the input arguments before the series is evaluated. The
+        map depends on the ``domain`` and ``window``; if the current
+        ``domain`` is equal to the ``window`` the resulting map is the
+        identity.  If the coeffients of the ``$name`` instance are to be
+        used by themselves outside this class, then the linear function
+        must be substituted for the ``x`` in the standard representation of
+        the base polynomials.
+
+        Returns
+        -------
+        off, scl : floats or complex
+            The mapping function is defined by ``off + scl*x``.
+
+        Notes
+        -----
+        If the current domain is the interval ``[l_1, r_1]`` and the window
+        is ``[l_2, r_2]``, then the linear mapping function ``L`` is
+        defined by the equations::
+
+            L(l_1) = l_2
+            L(r_1) = r_2
 
         """
-        return self.__class__(self.coef, self.domain)
+        return pu.mapparms(self.domain, self.window)
 
     def integ(self, m=1, k=[], lbnd=None) :
         """Integrate.
@@ -501,7 +608,7 @@ class $name(pu.PolyBase) :
         else :
             lbnd = off + scl*lbnd
         coef = ${nick}int(self.coef, m, k, lbnd, 1./scl)
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def deriv(self, m=1):
         """Differentiate.
@@ -527,7 +634,7 @@ class $name(pu.PolyBase) :
         """
         off, scl = self.mapparms()
         coef = ${nick}der(self.coef, m, scl)
-        return self.__class__(coef, self.domain)
+        return self.__class__(coef, self.domain, self.window)
 
     def roots(self) :
         """Return list of roots.
@@ -543,9 +650,9 @@ class $name(pu.PolyBase) :
 
         """
         roots = ${nick}roots(self.coef)
-        return pu.mapdomain(roots, $domain, self.domain)
+        return pu.mapdomain(roots, self.window, self.domain)
 
-    def linspace(self, n=100):
+    def linspace(self, n=100, domain=None):
         """Return x,y values at equally spaced points in domain.
 
         Returns x, y values at `n` equally spaced points across domain.
@@ -566,14 +673,17 @@ class $name(pu.PolyBase) :
         .. versionadded:: 1.5.0
 
         """
-        x = np.linspace(self.domain[0], self.domain[1], n)
+        if domain is None:
+            domain = self.domain
+        x = np.linspace(domain[0], domain[1], n)
         y = self(x)
         return x, y
 
 
 
     @staticmethod
-    def fit(x, y, deg, domain=None, rcond=None, full=False, w=None) :
+    def fit(x, y, deg, domain=None, rcond=None, full=False, w=None,
+        window=$domain):
         """Least squares fit to data.
 
         Return a `$name` instance that is the least squares fit to the data
@@ -616,6 +726,10 @@ class $name(pu.PolyBase) :
             ``w[i]*y[i]`` all have the same variance.  The default value is
             None.
             .. versionadded:: 1.5.0
+        window : {[beg, end]}, optional
+            Window to use for the returned $name instance. The default
+            value is ``$domain``
+            .. versionadded:: 1.6.0
 
         Returns
         -------
@@ -634,21 +748,25 @@ class $name(pu.PolyBase) :
         ${nick}fit : similar function
 
         """
-        if domain is None :
+        if domain is None:
             domain = pu.getdomain(x)
-        elif domain == [] :
+        elif domain == []:
             domain = $domain
-        xnew = pu.mapdomain(x, domain, $domain)
+
+        if window == []:
+            window = $domain
+
+        xnew = pu.mapdomain(x, domain, window)
         res = ${nick}fit(xnew, y, deg, w=w, rcond=rcond, full=full)
         if full :
             [coef, status] = res
-            return $name(coef, domain=domain), status
+            return $name(coef, domain=domain, window=window), status
         else :
             coef = res
-            return $name(coef, domain=domain)
+            return $name(coef, domain=domain, window=window)
 
     @staticmethod
-    def fromroots(roots, domain=$domain) :
+    def fromroots(roots, domain=$domain, window=$domain) :
         """Return $name object with specified roots.
 
         See ${nick}fromroots for full documentation.
@@ -660,12 +778,12 @@ class $name(pu.PolyBase) :
         """
         if domain is None :
             domain = pu.getdomain(roots)
-        rnew = pu.mapdomain(roots, domain, $domain)
+        rnew = pu.mapdomain(roots, domain, window)
         coef = ${nick}fromroots(rnew)
-        return $name(coef, domain=domain)
+        return $name(coef, domain=domain, window=window)
 
     @staticmethod
-    def identity(domain=$domain) :
+    def identity(domain=$domain, window=$domain) :
         """Identity function.
 
         If ``p`` is the returned $name object, then ``p(x) == x`` for all
@@ -676,13 +794,16 @@ class $name(pu.PolyBase) :
         domain : array_like
             The resulting array must be if the form ``[beg, end]``, where
             ``beg`` and ``end`` are the endpoints of the domain.
+        window : array_like
+            The resulting array must be if the form ``[beg, end]``, where
+            ``beg`` and ``end`` are the endpoints of the window.
 
         Returns
         -------
         identity : $name object
 
         """
-        off, scl = pu.mapparms($domain, domain)
+        off, scl = pu.mapparms(window, domain)
         coef = ${nick}line(off, scl)
-        return $name(coef, domain)
+        return $name(coef, domain, window)
 '''.replace('REL_IMPORT', rel_import))
diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py
index 65bb877f4..21e4728bf 100644
--- a/numpy/polynomial/tests/test_chebyshev.py
+++ b/numpy/polynomial/tests/test_chebyshev.py
@@ -462,7 +462,7 @@ class TestChebyshevClass(TestCase) :
     def test_degree(self) :
         assert_equal(self.p1.degree(), 2)
 
-    def test_trimdeg(self) :
+    def test_cutdeg(self) :
         assert_raises(ValueError, self.p1.cutdeg, .5)
         assert_raises(ValueError, self.p1.cutdeg, -1)
         assert_equal(len(self.p1.cutdeg(3)), 3)
@@ -564,3 +564,7 @@ class TestChebyshevClass(TestCase) :
         assert_almost_equal(p(x), x)
         p = ch.Chebyshev.identity([1,3])
         assert_almost_equal(p(x), x)
+#
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/numpy/polynomial/tests/test_hermite.py b/numpy/polynomial/tests/test_hermite.py
new file mode 100644
index 000000000..dea32f24a
--- /dev/null
+++ b/numpy/polynomial/tests/test_hermite.py
@@ -0,0 +1,537 @@
+"""Tests for hermendre module.
+
+"""
+from __future__ import division
+
+import numpy as np
+import numpy.polynomial.hermite as herm
+import numpy.polynomial.polynomial as poly
+from numpy.testing import *
+
+H0 = np.array([   1])
+H1 = np.array([0,     2])
+H2 = np.array([  -2, 0,      4])
+H3 = np.array([0,   -12, 0,      8])
+H4 = np.array([  12, 0,    -48, 0,    16])
+H5 = np.array([0,   120, 0,   -160, 0,    32])
+H6 = np.array([-120, 0,    720, 0,  -480, 0,    64])
+H7 = np.array([0, -1680, 0,   3360, 0, -1344, 0,   128])
+H8 = np.array([1680, 0, -13440, 0, 13440, 0, -3584, 0, 256])
+H9 = np.array([0, 30240, 0, -80640, 0, 48384, 0, -9216, 0, 512])
+
+Hlist = [H0, H1, H2, H3, H4, H5, H6, H7, H8, H9]
+
+
+def trim(x) :
+    return herm.hermtrim(x, tol=1e-6)
+
+
+class TestConstants(TestCase) :
+
+    def test_hermdomain(self) :
+        assert_equal(herm.hermdomain, [-1, 1])
+
+    def test_hermzero(self) :
+        assert_equal(herm.hermzero, [0])
+
+    def test_hermone(self) :
+        assert_equal(herm.hermone, [1])
+
+    def test_hermx(self) :
+        assert_equal(herm.hermx, [0, .5])
+
+
+class TestArithmetic(TestCase) :
+    x = np.linspace(-3, 3, 100)
+    y0 = poly.polyval(x, H0)
+    y1 = poly.polyval(x, H1)
+    y2 = poly.polyval(x, H2)
+    y3 = poly.polyval(x, H3)
+    y4 = poly.polyval(x, H4)
+    y5 = poly.polyval(x, H5)
+    y6 = poly.polyval(x, H6)
+    y7 = poly.polyval(x, H7)
+    y8 = poly.polyval(x, H8)
+    y9 = poly.polyval(x, H9)
+    y = [y0, y1, y2, y3, y4, y5, y6, y7, y8, y9]
+
+    def test_hermval(self) :
+        def f(x) :
+            return x*(x**2 - 1)
+
+        #check empty input
+        assert_equal(herm.hermval([], [1]).size, 0)
+
+        #check normal input)
+        for i in range(10) :
+            msg = "At i=%d" % i
+            ser = np.zeros
+            tgt = self.y[i]
+            res = herm.hermval(self.x, [0]*i + [1])
+            assert_almost_equal(res, tgt, err_msg=msg)
+
+        #check that shape is preserved
+        for i in range(3) :
+            dims = [2]*i
+            x = np.zeros(dims)
+            assert_equal(herm.hermval(x, [1]).shape, dims)
+            assert_equal(herm.hermval(x, [1,0]).shape, dims)
+            assert_equal(herm.hermval(x, [1,0,0]).shape, dims)
+
+    def test_hermadd(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                tgt = np.zeros(max(i,j) + 1)
+                tgt[i] += 1
+                tgt[j] += 1
+                res = herm.hermadd([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_hermsub(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                tgt = np.zeros(max(i,j) + 1)
+                tgt[i] += 1
+                tgt[j] -= 1
+                res = herm.hermsub([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_hermmulx(self):
+        assert_equal(herm.hermmulx([0]), [0])
+        assert_equal(herm.hermmulx([1]), [0,.5])
+        for i in range(1, 5):
+            ser = [0]*i + [1]
+            tgt = [0]*(i - 1) + [i, 0, .5]
+            assert_equal(herm.hermmulx(ser), tgt)
+
+    def test_hermmul(self) :
+        # check values of result
+        for i in range(5) :
+            pol1 = [0]*i + [1]
+            val1 = herm.hermval(self.x, pol1)
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                pol2 = [0]*j + [1]
+                val2 = herm.hermval(self.x, pol2)
+                pol3 = herm.hermmul(pol1, pol2)
+                val3 = herm.hermval(self.x, pol3)
+                assert_(len(pol3) == i + j + 1, msg)
+                assert_almost_equal(val3, val1*val2, err_msg=msg)
+
+    def test_hermdiv(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                ci = [0]*i + [1]
+                cj = [0]*j + [1]
+                tgt = herm.hermadd(ci, cj)
+                quo, rem = herm.hermdiv(tgt, ci)
+                res = herm.hermadd(herm.hermmul(quo, ci), rem)
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+
+class TestCalculus(TestCase) :
+
+    def test_hermint(self) :
+        # check exceptions
+        assert_raises(ValueError, herm.hermint, [0], .5)
+        assert_raises(ValueError, herm.hermint, [0], -1)
+        assert_raises(ValueError, herm.hermint, [0], 1, [0,0])
+
+        # test integration of zero polynomial
+        for i in range(2, 5):
+            k = [0]*(i - 2) + [1]
+            res = herm.hermint([0], m=i, k=k)
+            assert_almost_equal(res, [0, .5])
+
+        # check single integration with integration constant
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            tgt = [i] + [0]*i + [1/scl]
+            hermpol = herm.poly2herm(pol)
+            hermint = herm.hermint(hermpol, m=1, k=[i])
+            res = herm.herm2poly(hermint)
+            assert_almost_equal(trim(res), trim(tgt))
+
+        # check single integration with integration constant and lbnd
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            hermpol = herm.poly2herm(pol)
+            hermint = herm.hermint(hermpol, m=1, k=[i], lbnd=-1)
+            assert_almost_equal(herm.hermval(-1, hermint), i)
+
+        # check single integration with integration constant and scaling
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            tgt = [i] + [0]*i + [2/scl]
+            hermpol = herm.poly2herm(pol)
+            hermint = herm.hermint(hermpol, m=1, k=[i], scl=2)
+            res = herm.herm2poly(hermint)
+            assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with default k
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herm.hermint(tgt, m=1)
+                res = herm.hermint(pol, m=j)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with defined k
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herm.hermint(tgt, m=1, k=[k])
+                res = herm.hermint(pol, m=j, k=range(j))
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with lbnd
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herm.hermint(tgt, m=1, k=[k], lbnd=-1)
+                res = herm.hermint(pol, m=j, k=range(j), lbnd=-1)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with scaling
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herm.hermint(tgt, m=1, k=[k], scl=2)
+                res = herm.hermint(pol, m=j, k=range(j), scl=2)
+                assert_almost_equal(trim(res), trim(tgt))
+
+    def test_hermder(self) :
+        # check exceptions
+        assert_raises(ValueError, herm.hermder, [0], .5)
+        assert_raises(ValueError, herm.hermder, [0], -1)
+
+        # check that zeroth deriviative does nothing
+        for i in range(5) :
+            tgt = [1] + [0]*i
+            res = herm.hermder(tgt, m=0)
+            assert_equal(trim(res), trim(tgt))
+
+        # check that derivation is the inverse of integration
+        for i in range(5) :
+            for j in range(2,5) :
+                tgt = [1] + [0]*i
+                res = herm.hermder(herm.hermint(tgt, m=j), m=j)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check derivation with scaling
+        for i in range(5) :
+            for j in range(2,5) :
+                tgt = [1] + [0]*i
+                res = herm.hermder(herm.hermint(tgt, m=j, scl=2), m=j, scl=.5)
+                assert_almost_equal(trim(res), trim(tgt))
+
+
+class TestMisc(TestCase) :
+
+    def test_hermfromroots(self) :
+        res = herm.hermfromroots([])
+        assert_almost_equal(trim(res), [1])
+        for i in range(1,5) :
+            roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2])
+            pol = herm.hermfromroots(roots)
+            res = herm.hermval(roots, pol)
+            tgt = 0
+            assert_(len(pol) == i + 1)
+            assert_almost_equal(herm.herm2poly(pol)[-1], 1)
+            assert_almost_equal(res, tgt)
+
+    def test_hermroots(self) :
+        assert_almost_equal(herm.hermroots([1]), [])
+        assert_almost_equal(herm.hermroots([1, 1]), [-.5])
+        for i in range(2,5) :
+            tgt = np.linspace(-1, 1, i)
+            res = herm.hermroots(herm.hermfromroots(tgt))
+            assert_almost_equal(trim(res), trim(tgt))
+
+    def test_hermvander(self) :
+        # check for 1d x
+        x = np.arange(3)
+        v = herm.hermvander(x, 3)
+        assert_(v.shape == (3,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], herm.hermval(x, coef))
+
+        # check for 2d x
+        x = np.array([[1,2],[3,4],[5,6]])
+        v = herm.hermvander(x, 3)
+        assert_(v.shape == (3,2,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], herm.hermval(x, coef))
+
+    def test_hermfit(self) :
+        def f(x) :
+            return x*(x - 1)*(x - 2)
+
+        # Test exceptions
+        assert_raises(ValueError, herm.hermfit, [1],    [1],     -1)
+        assert_raises(TypeError,  herm.hermfit, [[1]],  [1],      0)
+        assert_raises(TypeError,  herm.hermfit, [],     [1],      0)
+        assert_raises(TypeError,  herm.hermfit, [1],    [[[1]]],  0)
+        assert_raises(TypeError,  herm.hermfit, [1, 2], [1],      0)
+        assert_raises(TypeError,  herm.hermfit, [1],    [1, 2],   0)
+        assert_raises(TypeError,  herm.hermfit, [1],    [1],   0, w=[[1]])
+        assert_raises(TypeError,  herm.hermfit, [1],    [1],   0, w=[1,1])
+
+        # Test fit
+        x = np.linspace(0,2)
+        y = f(x)
+        #
+        coef3 = herm.hermfit(x, y, 3)
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(herm.hermval(x, coef3), y)
+        #
+        coef4 = herm.hermfit(x, y, 4)
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(herm.hermval(x, coef4), y)
+        #
+        coef2d = herm.hermfit(x, np.array([y,y]).T, 3)
+        assert_almost_equal(coef2d, np.array([coef3,coef3]).T)
+        # test weighting
+        w = np.zeros_like(x)
+        yw = y.copy()
+        w[1::2] = 1
+        y[0::2] = 0
+        wcoef3 = herm.hermfit(x, yw, 3, w=w)
+        assert_almost_equal(wcoef3, coef3)
+        #
+        wcoef2d = herm.hermfit(x, np.array([yw,yw]).T, 3, w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3,coef3]).T)
+
+    def test_hermtrim(self) :
+        coef = [2, -1, 1, 0]
+
+        # Test exceptions
+        assert_raises(ValueError, herm.hermtrim, coef, -1)
+
+        # Test results
+        assert_equal(herm.hermtrim(coef), coef[:-1])
+        assert_equal(herm.hermtrim(coef, 1), coef[:-3])
+        assert_equal(herm.hermtrim(coef, 2), [0])
+
+    def test_hermline(self) :
+        assert_equal(herm.hermline(3,4), [3, 2])
+
+    def test_herm2poly(self) :
+        for i in range(10) :
+            assert_almost_equal(herm.herm2poly([0]*i + [1]), Hlist[i])
+
+    def test_poly2herm(self) :
+        for i in range(10) :
+            assert_almost_equal(herm.poly2herm(Hlist[i]), [0]*i + [1])
+
+
+def assert_poly_almost_equal(p1, p2):
+    assert_almost_equal(p1.coef, p2.coef)
+    assert_equal(p1.domain, p2.domain)
+
+
+class TestHermiteClass(TestCase) :
+
+    p1 = herm.Hermite([1,2,3])
+    p2 = herm.Hermite([1,2,3], [0,1])
+    p3 = herm.Hermite([1,2])
+    p4 = herm.Hermite([2,2,3])
+    p5 = herm.Hermite([3,2,3])
+
+    def test_equal(self) :
+        assert_(self.p1 == self.p1)
+        assert_(self.p2 == self.p2)
+        assert_(not self.p1 == self.p2)
+        assert_(not self.p1 == self.p3)
+        assert_(not self.p1 == [1,2,3])
+
+    def test_not_equal(self) :
+        assert_(not self.p1 != self.p1)
+        assert_(not self.p2 != self.p2)
+        assert_(self.p1 != self.p2)
+        assert_(self.p1 != self.p3)
+        assert_(self.p1 != [1,2,3])
+
+    def test_add(self) :
+        tgt = herm.Hermite([2,4,6])
+        assert_(self.p1 + self.p1 == tgt)
+        assert_(self.p1 + [1,2,3] == tgt)
+        assert_([1,2,3] + self.p1 == tgt)
+
+    def test_sub(self) :
+        tgt = herm.Hermite([1])
+        assert_(self.p4 - self.p1 == tgt)
+        assert_(self.p4 - [1,2,3] == tgt)
+        assert_([2,2,3] - self.p1 == tgt)
+
+    def test_mul(self) :
+        tgt = herm.Hermite([ 81.,  52.,  82.,  12.,   9.])
+        assert_poly_almost_equal(self.p1 * self.p1, tgt)
+        assert_poly_almost_equal(self.p1 * [1,2,3], tgt)
+        assert_poly_almost_equal([1,2,3] * self.p1, tgt)
+
+    def test_floordiv(self) :
+        tgt = herm.Hermite([1])
+        assert_(self.p4 // self.p1 == tgt)
+        assert_(self.p4 // [1,2,3] == tgt)
+        assert_([2,2,3] // self.p1 == tgt)
+
+    def test_mod(self) :
+        tgt = herm.Hermite([1])
+        assert_((self.p4 % self.p1) == tgt)
+        assert_((self.p4 % [1,2,3]) == tgt)
+        assert_(([2,2,3] % self.p1) == tgt)
+
+    def test_divmod(self) :
+        tquo = herm.Hermite([1])
+        trem = herm.Hermite([2])
+        quo, rem = divmod(self.p5, self.p1)
+        assert_(quo == tquo and rem == trem)
+        quo, rem = divmod(self.p5, [1,2,3])
+        assert_(quo == tquo and rem == trem)
+        quo, rem = divmod([3,2,3], self.p1)
+        assert_(quo == tquo and rem == trem)
+
+    def test_pow(self) :
+        tgt = herm.Hermite([1])
+        for i in range(5) :
+            res = self.p1**i
+            assert_(res == tgt)
+            tgt = tgt*self.p1
+
+    def test_call(self) :
+        # domain = [-1, 1]
+        x = np.linspace(-1, 1)
+        tgt = 3*(4*x**2 - 2) + 2*(2*x) + 1
+        assert_almost_equal(self.p1(x), tgt)
+
+        # domain = [0, 1]
+        x = np.linspace(0, 1)
+        xx = 2*x - 1
+        assert_almost_equal(self.p2(x), self.p1(xx))
+
+    def test_degree(self) :
+        assert_equal(self.p1.degree(), 2)
+
+    def test_cutdeg(self) :
+        assert_raises(ValueError, self.p1.cutdeg, .5)
+        assert_raises(ValueError, self.p1.cutdeg, -1)
+        assert_equal(len(self.p1.cutdeg(3)), 3)
+        assert_equal(len(self.p1.cutdeg(2)), 3)
+        assert_equal(len(self.p1.cutdeg(1)), 2)
+        assert_equal(len(self.p1.cutdeg(0)), 1)
+
+    def test_convert(self) :
+        x = np.linspace(-1,1)
+        p = self.p1.convert(domain=[0,1])
+        assert_almost_equal(p(x), self.p1(x))
+
+    def test_mapparms(self) :
+        parms = self.p2.mapparms()
+        assert_almost_equal(parms, [-1, 2])
+
+    def test_trim(self) :
+        coef = [1, 1e-6, 1e-12, 0]
+        p = herm.Hermite(coef)
+        assert_equal(p.trim().coef, coef[:3])
+        assert_equal(p.trim(1e-10).coef, coef[:2])
+        assert_equal(p.trim(1e-5).coef, coef[:1])
+
+    def test_truncate(self) :
+        assert_raises(ValueError, self.p1.truncate, .5)
+        assert_raises(ValueError, self.p1.truncate, 0)
+        assert_equal(len(self.p1.truncate(4)), 3)
+        assert_equal(len(self.p1.truncate(3)), 3)
+        assert_equal(len(self.p1.truncate(2)), 2)
+        assert_equal(len(self.p1.truncate(1)), 1)
+
+    def test_copy(self) :
+        p = self.p1.copy()
+        assert_(self.p1 == p)
+
+    def test_integ(self) :
+        p = self.p2.integ()
+        assert_almost_equal(p.coef, herm.hermint([1,2,3], 1, 0, scl=.5))
+        p = self.p2.integ(lbnd=0)
+        assert_almost_equal(p(0), 0)
+        p = self.p2.integ(1, 1)
+        assert_almost_equal(p.coef, herm.hermint([1,2,3], 1, 1, scl=.5))
+        p = self.p2.integ(2, [1, 2])
+        assert_almost_equal(p.coef, herm.hermint([1,2,3], 2, [1,2], scl=.5))
+
+    def test_deriv(self) :
+        p = self.p2.integ(2, [1, 2])
+        assert_almost_equal(p.deriv(1).coef, self.p2.integ(1, [1]).coef)
+        assert_almost_equal(p.deriv(2).coef, self.p2.coef)
+
+    def test_roots(self) :
+        p = herm.Hermite(herm.poly2herm([0, -1, 0, 1]), [0, 1])
+        res = p.roots()
+        tgt = [0, .5, 1]
+        assert_almost_equal(res, tgt)
+
+    def test_linspace(self):
+        xdes = np.linspace(0, 1, 20)
+        ydes = self.p2(xdes)
+        xres, yres = self.p2.linspace(20)
+        assert_almost_equal(xres, xdes)
+        assert_almost_equal(yres, ydes)
+
+    def test_fromroots(self) :
+        roots = [0, .5, 1]
+        p = herm.Hermite.fromroots(roots, domain=[0, 1])
+        res = p.coef
+        tgt = herm.poly2herm([0, -1, 0, 1])
+        assert_almost_equal(res, tgt)
+
+    def test_fit(self) :
+        def f(x) :
+            return x*(x - 1)*(x - 2)
+        x = np.linspace(0,3)
+        y = f(x)
+
+        # test default value of domain
+        p = herm.Hermite.fit(x, y, 3)
+        assert_almost_equal(p.domain, [0,3])
+
+        # test that fit works in given domains
+        p = herm.Hermite.fit(x, y, 3, None)
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [0,3])
+        p = herm.Hermite.fit(x, y, 3, [])
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [-1, 1])
+        # test that fit accepts weights.
+        w = np.zeros_like(x)
+        yw = y.copy()
+        w[1::2] = 1
+        yw[0::2] = 0
+        p = herm.Hermite.fit(x, yw, 3, w=w)
+        assert_almost_equal(p(x), y)
+
+    def test_identity(self) :
+        x = np.linspace(0,3)
+        p = herm.Hermite.identity()
+        assert_almost_equal(p(x), x)
+        p = herm.Hermite.identity([1,3])
+        assert_almost_equal(p(x), x)
+#
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py
new file mode 100644
index 000000000..6026e0e64
--- /dev/null
+++ b/numpy/polynomial/tests/test_hermite_e.py
@@ -0,0 +1,536 @@
+"""Tests for hermeendre module.
+
+"""
+from __future__ import division
+
+import numpy as np
+import numpy.polynomial.hermite_e as herme
+import numpy.polynomial.polynomial as poly
+from numpy.testing import *
+
+He0 = np.array([ 1 ])
+He1 = np.array([ 0 , 1 ])
+He2 = np.array([ -1 ,0 , 1 ])
+He3 = np.array([ 0 , -3 ,0 , 1 ])
+He4 = np.array([ 3 ,0 , -6 ,0 , 1 ])
+He5 = np.array([ 0 , 15 ,0 , -10 ,0 , 1 ])
+He6 = np.array([ -15 ,0 , 45 ,0 , -15 ,0 , 1 ])
+He7 = np.array([ 0 , -105 ,0 , 105 ,0 , -21 ,0 , 1 ])
+He8 = np.array([ 105 ,0 , -420 ,0 , 210 ,0 , -28 ,0 , 1 ])
+He9 = np.array([ 0 , 945 ,0 , -1260 ,0 , 378 ,0 , -36 ,0 , 1 ])
+
+Helist = [He0, He1, He2, He3, He4, He5, He6, He7, He8, He9]
+
+def trim(x) :
+    return herme.hermetrim(x, tol=1e-6)
+
+
+class TestConstants(TestCase) :
+
+    def test_hermedomain(self) :
+        assert_equal(herme.hermedomain, [-1, 1])
+
+    def test_hermezero(self) :
+        assert_equal(herme.hermezero, [0])
+
+    def test_hermeone(self) :
+        assert_equal(herme.hermeone, [1])
+
+    def test_hermex(self) :
+        assert_equal(herme.hermex, [0, 1])
+
+
+class TestArithmetic(TestCase) :
+    x = np.linspace(-3, 3, 100)
+    y0 = poly.polyval(x, He0)
+    y1 = poly.polyval(x, He1)
+    y2 = poly.polyval(x, He2)
+    y3 = poly.polyval(x, He3)
+    y4 = poly.polyval(x, He4)
+    y5 = poly.polyval(x, He5)
+    y6 = poly.polyval(x, He6)
+    y7 = poly.polyval(x, He7)
+    y8 = poly.polyval(x, He8)
+    y9 = poly.polyval(x, He9)
+    y = [y0, y1, y2, y3, y4, y5, y6, y7, y8, y9]
+
+    def test_hermeval(self) :
+        def f(x) :
+            return x*(x**2 - 1)
+
+        #check empty input
+        assert_equal(herme.hermeval([], [1]).size, 0)
+
+        #check normal input)
+        for i in range(10) :
+            msg = "At i=%d" % i
+            ser = np.zeros
+            tgt = self.y[i]
+            res = herme.hermeval(self.x, [0]*i + [1])
+            assert_almost_equal(res, tgt, err_msg=msg)
+
+        #check that shape is preserved
+        for i in range(3) :
+            dims = [2]*i
+            x = np.zeros(dims)
+            assert_equal(herme.hermeval(x, [1]).shape, dims)
+            assert_equal(herme.hermeval(x, [1,0]).shape, dims)
+            assert_equal(herme.hermeval(x, [1,0,0]).shape, dims)
+
+    def test_hermeadd(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                tgt = np.zeros(max(i,j) + 1)
+                tgt[i] += 1
+                tgt[j] += 1
+                res = herme.hermeadd([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_hermesub(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                tgt = np.zeros(max(i,j) + 1)
+                tgt[i] += 1
+                tgt[j] -= 1
+                res = herme.hermesub([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_hermemulx(self):
+        assert_equal(herme.hermemulx([0]), [0])
+        assert_equal(herme.hermemulx([1]), [0,1])
+        for i in range(1, 5):
+            ser = [0]*i + [1]
+            tgt = [0]*(i - 1) + [i, 0, 1]
+            assert_equal(herme.hermemulx(ser), tgt)
+
+    def test_hermemul(self) :
+        # check values of result
+        for i in range(5) :
+            pol1 = [0]*i + [1]
+            val1 = herme.hermeval(self.x, pol1)
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                pol2 = [0]*j + [1]
+                val2 = herme.hermeval(self.x, pol2)
+                pol3 = herme.hermemul(pol1, pol2)
+                val3 = herme.hermeval(self.x, pol3)
+                assert_(len(pol3) == i + j + 1, msg)
+                assert_almost_equal(val3, val1*val2, err_msg=msg)
+
+    def test_hermediv(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                ci = [0]*i + [1]
+                cj = [0]*j + [1]
+                tgt = herme.hermeadd(ci, cj)
+                quo, rem = herme.hermediv(tgt, ci)
+                res = herme.hermeadd(herme.hermemul(quo, ci), rem)
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+
+class TestCalculus(TestCase) :
+
+    def test_hermeint(self) :
+        # check exceptions
+        assert_raises(ValueError, herme.hermeint, [0], .5)
+        assert_raises(ValueError, herme.hermeint, [0], -1)
+        assert_raises(ValueError, herme.hermeint, [0], 1, [0,0])
+
+        # test integration of zero polynomial
+        for i in range(2, 5):
+            k = [0]*(i - 2) + [1]
+            res = herme.hermeint([0], m=i, k=k)
+            assert_almost_equal(res, [0, 1])
+
+        # check single integration with integration constant
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            tgt = [i] + [0]*i + [1/scl]
+            hermepol = herme.poly2herme(pol)
+            hermeint = herme.hermeint(hermepol, m=1, k=[i])
+            res = herme.herme2poly(hermeint)
+            assert_almost_equal(trim(res), trim(tgt))
+
+        # check single integration with integration constant and lbnd
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            hermepol = herme.poly2herme(pol)
+            hermeint = herme.hermeint(hermepol, m=1, k=[i], lbnd=-1)
+            assert_almost_equal(herme.hermeval(-1, hermeint), i)
+
+        # check single integration with integration constant and scaling
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            tgt = [i] + [0]*i + [2/scl]
+            hermepol = herme.poly2herme(pol)
+            hermeint = herme.hermeint(hermepol, m=1, k=[i], scl=2)
+            res = herme.herme2poly(hermeint)
+            assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with default k
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herme.hermeint(tgt, m=1)
+                res = herme.hermeint(pol, m=j)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with defined k
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herme.hermeint(tgt, m=1, k=[k])
+                res = herme.hermeint(pol, m=j, k=range(j))
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with lbnd
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herme.hermeint(tgt, m=1, k=[k], lbnd=-1)
+                res = herme.hermeint(pol, m=j, k=range(j), lbnd=-1)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with scaling
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = herme.hermeint(tgt, m=1, k=[k], scl=2)
+                res = herme.hermeint(pol, m=j, k=range(j), scl=2)
+                assert_almost_equal(trim(res), trim(tgt))
+
+    def test_hermeder(self) :
+        # check exceptions
+        assert_raises(ValueError, herme.hermeder, [0], .5)
+        assert_raises(ValueError, herme.hermeder, [0], -1)
+
+        # check that zeroth deriviative does nothing
+        for i in range(5) :
+            tgt = [1] + [0]*i
+            res = herme.hermeder(tgt, m=0)
+            assert_equal(trim(res), trim(tgt))
+
+        # check that derivation is the inverse of integration
+        for i in range(5) :
+            for j in range(2,5) :
+                tgt = [1] + [0]*i
+                res = herme.hermeder(herme.hermeint(tgt, m=j), m=j)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check derivation with scaling
+        for i in range(5) :
+            for j in range(2,5) :
+                tgt = [1] + [0]*i
+                res = herme.hermeder(herme.hermeint(tgt, m=j, scl=2), m=j, scl=.5)
+                assert_almost_equal(trim(res), trim(tgt))
+
+
+class TestMisc(TestCase) :
+
+    def test_hermefromroots(self) :
+        res = herme.hermefromroots([])
+        assert_almost_equal(trim(res), [1])
+        for i in range(1,5) :
+            roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2])
+            pol = herme.hermefromroots(roots)
+            res = herme.hermeval(roots, pol)
+            tgt = 0
+            assert_(len(pol) == i + 1)
+            assert_almost_equal(herme.herme2poly(pol)[-1], 1)
+            assert_almost_equal(res, tgt)
+
+    def test_hermeroots(self) :
+        assert_almost_equal(herme.hermeroots([1]), [])
+        assert_almost_equal(herme.hermeroots([1, 1]), [-.5])
+        for i in range(2,5) :
+            tgt = np.linspace(-1, 1, i)
+            res = herme.hermeroots(herme.hermefromroots(tgt))
+            assert_almost_equal(trim(res), trim(tgt))
+
+    def test_hermevander(self) :
+        # check for 1d x
+        x = np.arange(3)
+        v = herme.hermevander(x, 3)
+        assert_(v.shape == (3,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], herme.hermeval(x, coef))
+
+        # check for 2d x
+        x = np.array([[1,2],[3,4],[5,6]])
+        v = herme.hermevander(x, 3)
+        assert_(v.shape == (3,2,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], herme.hermeval(x, coef))
+
+    def test_hermefit(self) :
+        def f(x) :
+            return x*(x - 1)*(x - 2)
+
+        # Test exceptions
+        assert_raises(ValueError, herme.hermefit, [1],    [1],     -1)
+        assert_raises(TypeError,  herme.hermefit, [[1]],  [1],      0)
+        assert_raises(TypeError,  herme.hermefit, [],     [1],      0)
+        assert_raises(TypeError,  herme.hermefit, [1],    [[[1]]],  0)
+        assert_raises(TypeError,  herme.hermefit, [1, 2], [1],      0)
+        assert_raises(TypeError,  herme.hermefit, [1],    [1, 2],   0)
+        assert_raises(TypeError,  herme.hermefit, [1],    [1],   0, w=[[1]])
+        assert_raises(TypeError,  herme.hermefit, [1],    [1],   0, w=[1,1])
+
+        # Test fit
+        x = np.linspace(0,2)
+        y = f(x)
+        #
+        coef3 = herme.hermefit(x, y, 3)
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(herme.hermeval(x, coef3), y)
+        #
+        coef4 = herme.hermefit(x, y, 4)
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(herme.hermeval(x, coef4), y)
+        #
+        coef2d = herme.hermefit(x, np.array([y,y]).T, 3)
+        assert_almost_equal(coef2d, np.array([coef3,coef3]).T)
+        # test weighting
+        w = np.zeros_like(x)
+        yw = y.copy()
+        w[1::2] = 1
+        y[0::2] = 0
+        wcoef3 = herme.hermefit(x, yw, 3, w=w)
+        assert_almost_equal(wcoef3, coef3)
+        #
+        wcoef2d = herme.hermefit(x, np.array([yw,yw]).T, 3, w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3,coef3]).T)
+
+    def test_hermetrim(self) :
+        coef = [2, -1, 1, 0]
+
+        # Test exceptions
+        assert_raises(ValueError, herme.hermetrim, coef, -1)
+
+        # Test results
+        assert_equal(herme.hermetrim(coef), coef[:-1])
+        assert_equal(herme.hermetrim(coef, 1), coef[:-3])
+        assert_equal(herme.hermetrim(coef, 2), [0])
+
+    def test_hermeline(self) :
+        assert_equal(herme.hermeline(3,4), [3, 4])
+
+    def test_herme2poly(self) :
+        for i in range(10) :
+            assert_almost_equal(herme.herme2poly([0]*i + [1]), Helist[i])
+
+    def test_poly2herme(self) :
+        for i in range(10) :
+            assert_almost_equal(herme.poly2herme(Helist[i]), [0]*i + [1])
+
+
+def assert_poly_almost_equal(p1, p2):
+    assert_almost_equal(p1.coef, p2.coef)
+    assert_equal(p1.domain, p2.domain)
+
+
+class TestHermiteEClass(TestCase) :
+
+    p1 = herme.HermiteE([1,2,3])
+    p2 = herme.HermiteE([1,2,3], [0,1])
+    p3 = herme.HermiteE([1,2])
+    p4 = herme.HermiteE([2,2,3])
+    p5 = herme.HermiteE([3,2,3])
+
+    def test_equal(self) :
+        assert_(self.p1 == self.p1)
+        assert_(self.p2 == self.p2)
+        assert_(not self.p1 == self.p2)
+        assert_(not self.p1 == self.p3)
+        assert_(not self.p1 == [1,2,3])
+
+    def test_not_equal(self) :
+        assert_(not self.p1 != self.p1)
+        assert_(not self.p2 != self.p2)
+        assert_(self.p1 != self.p2)
+        assert_(self.p1 != self.p3)
+        assert_(self.p1 != [1,2,3])
+
+    def test_add(self) :
+        tgt = herme.HermiteE([2,4,6])
+        assert_(self.p1 + self.p1 == tgt)
+        assert_(self.p1 + [1,2,3] == tgt)
+        assert_([1,2,3] + self.p1 == tgt)
+
+    def test_sub(self) :
+        tgt = herme.HermiteE([1])
+        assert_(self.p4 - self.p1 == tgt)
+        assert_(self.p4 - [1,2,3] == tgt)
+        assert_([2,2,3] - self.p1 == tgt)
+
+    def test_mul(self) :
+        tgt = herme.HermiteE([ 23.,  28.,  46.,  12.,   9.])
+        assert_poly_almost_equal(self.p1 * self.p1, tgt)
+        assert_poly_almost_equal(self.p1 * [1,2,3], tgt)
+        assert_poly_almost_equal([1,2,3] * self.p1, tgt)
+
+    def test_floordiv(self) :
+        tgt = herme.HermiteE([1])
+        assert_(self.p4 // self.p1 == tgt)
+        assert_(self.p4 // [1,2,3] == tgt)
+        assert_([2,2,3] // self.p1 == tgt)
+
+    def test_mod(self) :
+        tgt = herme.HermiteE([1])
+        assert_((self.p4 % self.p1) == tgt)
+        assert_((self.p4 % [1,2,3]) == tgt)
+        assert_(([2,2,3] % self.p1) == tgt)
+
+    def test_divmod(self) :
+        tquo = herme.HermiteE([1])
+        trem = herme.HermiteE([2])
+        quo, rem = divmod(self.p5, self.p1)
+        assert_(quo == tquo and rem == trem)
+        quo, rem = divmod(self.p5, [1,2,3])
+        assert_(quo == tquo and rem == trem)
+        quo, rem = divmod([3,2,3], self.p1)
+        assert_(quo == tquo and rem == trem)
+
+    def test_pow(self) :
+        tgt = herme.HermiteE([1])
+        for i in range(5) :
+            res = self.p1**i
+            assert_(res == tgt)
+            tgt = tgt*self.p1
+
+    def test_call(self) :
+        # domain = [-1, 1]
+        x = np.linspace(-1, 1)
+        tgt = 3*(x**2 - 1) + 2*(x) + 1
+        assert_almost_equal(self.p1(x), tgt)
+
+        # domain = [0, 1]
+        x = np.linspace(0, 1)
+        xx = 2*x - 1
+        assert_almost_equal(self.p2(x), self.p1(xx))
+
+    def test_degree(self) :
+        assert_equal(self.p1.degree(), 2)
+
+    def test_cutdeg(self) :
+        assert_raises(ValueError, self.p1.cutdeg, .5)
+        assert_raises(ValueError, self.p1.cutdeg, -1)
+        assert_equal(len(self.p1.cutdeg(3)), 3)
+        assert_equal(len(self.p1.cutdeg(2)), 3)
+        assert_equal(len(self.p1.cutdeg(1)), 2)
+        assert_equal(len(self.p1.cutdeg(0)), 1)
+
+    def test_convert(self) :
+        x = np.linspace(-1,1)
+        p = self.p1.convert(domain=[0,1])
+        assert_almost_equal(p(x), self.p1(x))
+
+    def test_mapparms(self) :
+        parms = self.p2.mapparms()
+        assert_almost_equal(parms, [-1, 2])
+
+    def test_trim(self) :
+        coef = [1, 1e-6, 1e-12, 0]
+        p = herme.HermiteE(coef)
+        assert_equal(p.trim().coef, coef[:3])
+        assert_equal(p.trim(1e-10).coef, coef[:2])
+        assert_equal(p.trim(1e-5).coef, coef[:1])
+
+    def test_truncate(self) :
+        assert_raises(ValueError, self.p1.truncate, .5)
+        assert_raises(ValueError, self.p1.truncate, 0)
+        assert_equal(len(self.p1.truncate(4)), 3)
+        assert_equal(len(self.p1.truncate(3)), 3)
+        assert_equal(len(self.p1.truncate(2)), 2)
+        assert_equal(len(self.p1.truncate(1)), 1)
+
+    def test_copy(self) :
+        p = self.p1.copy()
+        assert_(self.p1 == p)
+
+    def test_integ(self) :
+        p = self.p2.integ()
+        assert_almost_equal(p.coef, herme.hermeint([1,2,3], 1, 0, scl=.5))
+        p = self.p2.integ(lbnd=0)
+        assert_almost_equal(p(0), 0)
+        p = self.p2.integ(1, 1)
+        assert_almost_equal(p.coef, herme.hermeint([1,2,3], 1, 1, scl=.5))
+        p = self.p2.integ(2, [1, 2])
+        assert_almost_equal(p.coef, herme.hermeint([1,2,3], 2, [1,2], scl=.5))
+
+    def test_deriv(self) :
+        p = self.p2.integ(2, [1, 2])
+        assert_almost_equal(p.deriv(1).coef, self.p2.integ(1, [1]).coef)
+        assert_almost_equal(p.deriv(2).coef, self.p2.coef)
+
+    def test_roots(self) :
+        p = herme.HermiteE(herme.poly2herme([0, -1, 0, 1]), [0, 1])
+        res = p.roots()
+        tgt = [0, .5, 1]
+        assert_almost_equal(res, tgt)
+
+    def test_linspace(self):
+        xdes = np.linspace(0, 1, 20)
+        ydes = self.p2(xdes)
+        xres, yres = self.p2.linspace(20)
+        assert_almost_equal(xres, xdes)
+        assert_almost_equal(yres, ydes)
+
+    def test_fromroots(self) :
+        roots = [0, .5, 1]
+        p = herme.HermiteE.fromroots(roots, domain=[0, 1])
+        res = p.coef
+        tgt = herme.poly2herme([0, -1, 0, 1])
+        assert_almost_equal(res, tgt)
+
+    def test_fit(self) :
+        def f(x) :
+            return x*(x - 1)*(x - 2)
+        x = np.linspace(0,3)
+        y = f(x)
+
+        # test default value of domain
+        p = herme.HermiteE.fit(x, y, 3)
+        assert_almost_equal(p.domain, [0,3])
+
+        # test that fit works in given domains
+        p = herme.HermiteE.fit(x, y, 3, None)
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [0,3])
+        p = herme.HermiteE.fit(x, y, 3, [])
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [-1, 1])
+        # test that fit accepts weights.
+        w = np.zeros_like(x)
+        yw = y.copy()
+        w[1::2] = 1
+        yw[0::2] = 0
+        p = herme.HermiteE.fit(x, yw, 3, w=w)
+        assert_almost_equal(p(x), y)
+
+    def test_identity(self) :
+        x = np.linspace(0,3)
+        p = herme.HermiteE.identity()
+        assert_almost_equal(p(x), x)
+        p = herme.HermiteE.identity([1,3])
+        assert_almost_equal(p(x), x)
+
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/numpy/polynomial/tests/test_laguerre.py b/numpy/polynomial/tests/test_laguerre.py
new file mode 100644
index 000000000..f3a4a930b
--- /dev/null
+++ b/numpy/polynomial/tests/test_laguerre.py
@@ -0,0 +1,530 @@
+"""Tests for hermendre module.
+
+"""
+from __future__ import division
+
+import numpy as np
+import numpy.polynomial.laguerre as lag
+import numpy.polynomial.polynomial as poly
+from numpy.testing import *
+
+L0 = np.array([1 ])/1
+L1 = np.array([1 , -1 ])/1
+L2 = np.array([2 , -4 , 1 ])/2
+L3 = np.array([6 , -18 , 9 , -1 ])/6
+L4 = np.array([24 , -96 , 72 , -16 , 1 ])/24
+L5 = np.array([120 , -600 , 600 , -200 , 25 , -1 ])/120
+L6 = np.array([720 , -4320 , 5400 , -2400 , 450 , -36 , 1 ])/720
+
+Llist = [L0, L1, L2, L3, L4, L5, L6]
+
+def trim(x) :
+    return lag.lagtrim(x, tol=1e-6)
+
+
+class TestConstants(TestCase) :
+
+    def test_lagdomain(self) :
+        assert_equal(lag.lagdomain, [0, 1])
+
+    def test_lagzero(self) :
+        assert_equal(lag.lagzero, [0])
+
+    def test_lagone(self) :
+        assert_equal(lag.lagone, [1])
+
+    def test_lagx(self) :
+        assert_equal(lag.lagx, [1, -1])
+
+
+class TestArithmetic(TestCase) :
+    x = np.linspace(-3, 3, 100)
+    y0 = poly.polyval(x, L0)
+    y1 = poly.polyval(x, L1)
+    y2 = poly.polyval(x, L2)
+    y3 = poly.polyval(x, L3)
+    y4 = poly.polyval(x, L4)
+    y5 = poly.polyval(x, L5)
+    y6 = poly.polyval(x, L6)
+    y = [y0, y1, y2, y3, y4, y5, y6]
+
+    def test_lagval(self) :
+        def f(x) :
+            return x*(x**2 - 1)
+
+        #check empty input
+        assert_equal(lag.lagval([], [1]).size, 0)
+
+        #check normal input)
+        for i in range(7) :
+            msg = "At i=%d" % i
+            ser = np.zeros
+            tgt = self.y[i]
+            res = lag.lagval(self.x, [0]*i + [1])
+            assert_almost_equal(res, tgt, err_msg=msg)
+
+        #check that shape is preserved
+        for i in range(3) :
+            dims = [2]*i
+            x = np.zeros(dims)
+            assert_equal(lag.lagval(x, [1]).shape, dims)
+            assert_equal(lag.lagval(x, [1,0]).shape, dims)
+            assert_equal(lag.lagval(x, [1,0,0]).shape, dims)
+
+    def test_lagadd(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                tgt = np.zeros(max(i,j) + 1)
+                tgt[i] += 1
+                tgt[j] += 1
+                res = lag.lagadd([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_lagsub(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                tgt = np.zeros(max(i,j) + 1)
+                tgt[i] += 1
+                tgt[j] -= 1
+                res = lag.lagsub([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_lagmulx(self):
+        assert_equal(lag.lagmulx([0]), [0])
+        assert_equal(lag.lagmulx([1]), [1,-1])
+        for i in range(1, 5):
+            ser = [0]*i + [1]
+            tgt = [0]*(i - 1) + [-i, 2*i + 1, -(i + 1)]
+            assert_almost_equal(lag.lagmulx(ser), tgt)
+
+    def test_lagmul(self) :
+        # check values of result
+        for i in range(5) :
+            pol1 = [0]*i + [1]
+            val1 = lag.lagval(self.x, pol1)
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                pol2 = [0]*j + [1]
+                val2 = lag.lagval(self.x, pol2)
+                pol3 = lag.lagmul(pol1, pol2)
+                val3 = lag.lagval(self.x, pol3)
+                assert_(len(pol3) == i + j + 1, msg)
+                assert_almost_equal(val3, val1*val2, err_msg=msg)
+
+    def test_lagdiv(self) :
+        for i in range(5) :
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                ci = [0]*i + [1]
+                cj = [0]*j + [1]
+                tgt = lag.lagadd(ci, cj)
+                quo, rem = lag.lagdiv(tgt, ci)
+                res = lag.lagadd(lag.lagmul(quo, ci), rem)
+                assert_almost_equal(trim(res), trim(tgt), err_msg=msg)
+
+
+class TestCalculus(TestCase) :
+
+    def test_lagint(self) :
+        # check exceptions
+        assert_raises(ValueError, lag.lagint, [0], .5)
+        assert_raises(ValueError, lag.lagint, [0], -1)
+        assert_raises(ValueError, lag.lagint, [0], 1, [0,0])
+
+        # test integration of zero polynomial
+        for i in range(2, 5):
+            k = [0]*(i - 2) + [1]
+            res = lag.lagint([0], m=i, k=k)
+            assert_almost_equal(res, [1, -1])
+
+        # check single integration with integration constant
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            tgt = [i] + [0]*i + [1/scl]
+            lagpol = lag.poly2lag(pol)
+            lagint = lag.lagint(lagpol, m=1, k=[i])
+            res = lag.lag2poly(lagint)
+            assert_almost_equal(trim(res), trim(tgt))
+
+        # check single integration with integration constant and lbnd
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            lagpol = lag.poly2lag(pol)
+            lagint = lag.lagint(lagpol, m=1, k=[i], lbnd=-1)
+            assert_almost_equal(lag.lagval(-1, lagint), i)
+
+        # check single integration with integration constant and scaling
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            tgt = [i] + [0]*i + [2/scl]
+            lagpol = lag.poly2lag(pol)
+            lagint = lag.lagint(lagpol, m=1, k=[i], scl=2)
+            res = lag.lag2poly(lagint)
+            assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with default k
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = lag.lagint(tgt, m=1)
+                res = lag.lagint(pol, m=j)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with defined k
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = lag.lagint(tgt, m=1, k=[k])
+                res = lag.lagint(pol, m=j, k=range(j))
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with lbnd
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = lag.lagint(tgt, m=1, k=[k], lbnd=-1)
+                res = lag.lagint(pol, m=j, k=range(j), lbnd=-1)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check multiple integrations with scaling
+        for i in range(5) :
+            for j in range(2,5) :
+                pol = [0]*i + [1]
+                tgt = pol[:]
+                for k in range(j) :
+                    tgt = lag.lagint(tgt, m=1, k=[k], scl=2)
+                res = lag.lagint(pol, m=j, k=range(j), scl=2)
+                assert_almost_equal(trim(res), trim(tgt))
+
+    def test_lagder(self) :
+        # check exceptions
+        assert_raises(ValueError, lag.lagder, [0], .5)
+        assert_raises(ValueError, lag.lagder, [0], -1)
+
+        # check that zeroth deriviative does nothing
+        for i in range(5) :
+            tgt = [1] + [0]*i
+            res = lag.lagder(tgt, m=0)
+            assert_equal(trim(res), trim(tgt))
+
+        # check that derivation is the inverse of integration
+        for i in range(5) :
+            for j in range(2,5) :
+                tgt = [1] + [0]*i
+                res = lag.lagder(lag.lagint(tgt, m=j), m=j)
+                assert_almost_equal(trim(res), trim(tgt))
+
+        # check derivation with scaling
+        for i in range(5) :
+            for j in range(2,5) :
+                tgt = [1] + [0]*i
+                res = lag.lagder(lag.lagint(tgt, m=j, scl=2), m=j, scl=.5)
+                assert_almost_equal(trim(res), trim(tgt))
+
+
+class TestMisc(TestCase) :
+
+    def test_lagfromroots(self) :
+        res = lag.lagfromroots([])
+        assert_almost_equal(trim(res), [1])
+        for i in range(1,5) :
+            roots = np.cos(np.linspace(-np.pi, 0, 2*i + 1)[1::2])
+            pol = lag.lagfromroots(roots)
+            res = lag.lagval(roots, pol)
+            tgt = 0
+            assert_(len(pol) == i + 1)
+            assert_almost_equal(lag.lag2poly(pol)[-1], 1)
+            assert_almost_equal(res, tgt)
+
+    def test_lagroots(self) :
+        assert_almost_equal(lag.lagroots([1]), [])
+        assert_almost_equal(lag.lagroots([0, 1]), [1])
+        for i in range(2,5) :
+            tgt = np.linspace(0, 3, i)
+            res = lag.lagroots(lag.lagfromroots(tgt))
+            assert_almost_equal(trim(res), trim(tgt))
+
+    def test_lagvander(self) :
+        # check for 1d x
+        x = np.arange(3)
+        v = lag.lagvander(x, 3)
+        assert_(v.shape == (3,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], lag.lagval(x, coef))
+
+        # check for 2d x
+        x = np.array([[1,2],[3,4],[5,6]])
+        v = lag.lagvander(x, 3)
+        assert_(v.shape == (3,2,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], lag.lagval(x, coef))
+
+    def test_lagfit(self) :
+        def f(x) :
+            return x*(x - 1)*(x - 2)
+
+        # Test exceptions
+        assert_raises(ValueError, lag.lagfit, [1],    [1],     -1)
+        assert_raises(TypeError,  lag.lagfit, [[1]],  [1],      0)
+        assert_raises(TypeError,  lag.lagfit, [],     [1],      0)
+        assert_raises(TypeError,  lag.lagfit, [1],    [[[1]]],  0)
+        assert_raises(TypeError,  lag.lagfit, [1, 2], [1],      0)
+        assert_raises(TypeError,  lag.lagfit, [1],    [1, 2],   0)
+        assert_raises(TypeError,  lag.lagfit, [1],    [1],   0, w=[[1]])
+        assert_raises(TypeError,  lag.lagfit, [1],    [1],   0, w=[1,1])
+
+        # Test fit
+        x = np.linspace(0,2)
+        y = f(x)
+        #
+        coef3 = lag.lagfit(x, y, 3)
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(lag.lagval(x, coef3), y)
+        #
+        coef4 = lag.lagfit(x, y, 4)
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(lag.lagval(x, coef4), y)
+        #
+        coef2d = lag.lagfit(x, np.array([y,y]).T, 3)
+        assert_almost_equal(coef2d, np.array([coef3,coef3]).T)
+        # test weighting
+        w = np.zeros_like(x)
+        yw = y.copy()
+        w[1::2] = 1
+        y[0::2] = 0
+        wcoef3 = lag.lagfit(x, yw, 3, w=w)
+        assert_almost_equal(wcoef3, coef3)
+        #
+        wcoef2d = lag.lagfit(x, np.array([yw,yw]).T, 3, w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3,coef3]).T)
+
+    def test_lagtrim(self) :
+        coef = [2, -1, 1, 0]
+
+        # Test exceptions
+        assert_raises(ValueError, lag.lagtrim, coef, -1)
+
+        # Test results
+        assert_equal(lag.lagtrim(coef), coef[:-1])
+        assert_equal(lag.lagtrim(coef, 1), coef[:-3])
+        assert_equal(lag.lagtrim(coef, 2), [0])
+
+    def test_lagline(self) :
+        assert_equal(lag.lagline(3,4), [7, -4])
+
+    def test_lag2poly(self) :
+        for i in range(7) :
+            assert_almost_equal(lag.lag2poly([0]*i + [1]), Llist[i])
+
+    def test_poly2lag(self) :
+        for i in range(7) :
+            assert_almost_equal(lag.poly2lag(Llist[i]), [0]*i + [1])
+
+
+def assert_poly_almost_equal(p1, p2):
+    assert_almost_equal(p1.coef, p2.coef)
+    assert_equal(p1.domain, p2.domain)
+
+
+class TestLaguerreClass(TestCase) :
+
+    p1 = lag.Laguerre([1,2,3])
+    p2 = lag.Laguerre([1,2,3], [0,1])
+    p3 = lag.Laguerre([1,2])
+    p4 = lag.Laguerre([2,2,3])
+    p5 = lag.Laguerre([3,2,3])
+
+    def test_equal(self) :
+        assert_(self.p1 == self.p1)
+        assert_(self.p2 == self.p2)
+        assert_(not self.p1 == self.p2)
+        assert_(not self.p1 == self.p3)
+        assert_(not self.p1 == [1,2,3])
+
+    def test_not_equal(self) :
+        assert_(not self.p1 != self.p1)
+        assert_(not self.p2 != self.p2)
+        assert_(self.p1 != self.p2)
+        assert_(self.p1 != self.p3)
+        assert_(self.p1 != [1,2,3])
+
+    def test_add(self) :
+        tgt = lag.Laguerre([2,4,6])
+        assert_(self.p1 + self.p1 == tgt)
+        assert_(self.p1 + [1,2,3] == tgt)
+        assert_([1,2,3] + self.p1 == tgt)
+
+    def test_sub(self) :
+        tgt = lag.Laguerre([1])
+        assert_(self.p4 - self.p1 == tgt)
+        assert_(self.p4 - [1,2,3] == tgt)
+        assert_([2,2,3] - self.p1 == tgt)
+
+    def test_mul(self) :
+        tgt = lag.Laguerre([ 14., -16.,  56., -72.,  54.])
+        assert_poly_almost_equal(self.p1 * self.p1, tgt)
+        assert_poly_almost_equal(self.p1 * [1,2,3], tgt)
+        assert_poly_almost_equal([1,2,3] * self.p1, tgt)
+
+    def test_floordiv(self) :
+        tgt = lag.Laguerre([1])
+        assert_(self.p4 // self.p1 == tgt)
+        assert_(self.p4 // [1,2,3] == tgt)
+        assert_([2,2,3] // self.p1 == tgt)
+
+    def test_mod(self) :
+        tgt = lag.Laguerre([1])
+        assert_((self.p4 % self.p1) == tgt)
+        assert_((self.p4 % [1,2,3]) == tgt)
+        assert_(([2,2,3] % self.p1) == tgt)
+
+    def test_divmod(self) :
+        tquo = lag.Laguerre([1])
+        trem = lag.Laguerre([2])
+        quo, rem = divmod(self.p5, self.p1)
+        assert_(quo == tquo and rem == trem)
+        quo, rem = divmod(self.p5, [1,2,3])
+        assert_(quo == tquo and rem == trem)
+        quo, rem = divmod([3,2,3], self.p1)
+        assert_(quo == tquo and rem == trem)
+
+    def test_pow(self) :
+        tgt = lag.Laguerre([1])
+        for i in range(5) :
+            res = self.p1**i
+            assert_(res == tgt)
+            tgt = tgt*self.p1
+
+    def test_call(self) :
+        # domain = [0, 1]
+        x = np.linspace(0, 1)
+        tgt = 3*(.5*x**2 - 2*x + 1) + 2*(-x + 1) + 1
+        assert_almost_equal(self.p1(x), tgt)
+
+        # domain = [0, 1]
+        x = np.linspace(.5, 1)
+        xx = 2*x - 1
+        assert_almost_equal(self.p2(x), self.p1(xx))
+
+    def test_degree(self) :
+        assert_equal(self.p1.degree(), 2)
+
+    def test_cutdeg(self) :
+        assert_raises(ValueError, self.p1.cutdeg, .5)
+        assert_raises(ValueError, self.p1.cutdeg, -1)
+        assert_equal(len(self.p1.cutdeg(3)), 3)
+        assert_equal(len(self.p1.cutdeg(2)), 3)
+        assert_equal(len(self.p1.cutdeg(1)), 2)
+        assert_equal(len(self.p1.cutdeg(0)), 1)
+
+    def test_convert(self) :
+        x = np.linspace(-1,1)
+        p = self.p1.convert(domain=[0,1])
+        assert_almost_equal(p(x), self.p1(x))
+
+    def test_mapparms(self) :
+        parms = self.p2.mapparms()
+        assert_almost_equal(parms, [-1, 2])
+
+    def test_trim(self) :
+        coef = [1, 1e-6, 1e-12, 0]
+        p = lag.Laguerre(coef)
+        assert_equal(p.trim().coef, coef[:3])
+        assert_equal(p.trim(1e-10).coef, coef[:2])
+        assert_equal(p.trim(1e-5).coef, coef[:1])
+
+    def test_truncate(self) :
+        assert_raises(ValueError, self.p1.truncate, .5)
+        assert_raises(ValueError, self.p1.truncate, 0)
+        assert_equal(len(self.p1.truncate(4)), 3)
+        assert_equal(len(self.p1.truncate(3)), 3)
+        assert_equal(len(self.p1.truncate(2)), 2)
+        assert_equal(len(self.p1.truncate(1)), 1)
+
+    def test_copy(self) :
+        p = self.p1.copy()
+        assert_(self.p1 == p)
+
+    def test_integ(self) :
+        p = self.p2.integ()
+        assert_almost_equal(p.coef, lag.lagint([1,2,3], 1, 0, scl=.5))
+        p = self.p2.integ(lbnd=0)
+        assert_almost_equal(p(0), 0)
+        p = self.p2.integ(1, 1)
+        assert_almost_equal(p.coef, lag.lagint([1,2,3], 1, 1, scl=.5))
+        p = self.p2.integ(2, [1, 2])
+        assert_almost_equal(p.coef, lag.lagint([1,2,3], 2, [1,2], scl=.5))
+
+    def test_deriv(self) :
+        p = self.p2.integ(2, [1, 2])
+        assert_almost_equal(p.deriv(1).coef, self.p2.integ(1, [1]).coef)
+        assert_almost_equal(p.deriv(2).coef, self.p2.coef)
+
+    def test_roots(self) :
+        p = lag.Laguerre(lag.poly2lag([0, -1, 0, 1]), [0, 1])
+        res = p.roots()
+        tgt = [0, .5, 1]
+        assert_almost_equal(res, tgt)
+
+    def test_linspace(self):
+        xdes = np.linspace(0, 1, 20)
+        ydes = self.p2(xdes)
+        xres, yres = self.p2.linspace(20)
+        assert_almost_equal(xres, xdes)
+        assert_almost_equal(yres, ydes)
+
+    def test_fromroots(self) :
+        roots = [0, .5, 1]
+        p = lag.Laguerre.fromroots(roots, domain=[0, 1])
+        res = p.coef
+        tgt = lag.poly2lag([0, -1, 0, 1])
+        assert_almost_equal(res, tgt)
+
+    def test_fit(self) :
+        def f(x) :
+            return x*(x - 1)*(x - 2)
+        x = np.linspace(0,3)
+        y = f(x)
+
+        # test default value of domain
+        p = lag.Laguerre.fit(x, y, 3)
+        assert_almost_equal(p.domain, [0,3])
+
+        # test that fit works in given domains
+        p = lag.Laguerre.fit(x, y, 3, None)
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [0,3])
+        p = lag.Laguerre.fit(x, y, 3, [])
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [-1, 1])
+        # test that fit accepts weights.
+        w = np.zeros_like(x)
+        yw = y.copy()
+        w[1::2] = 1
+        yw[0::2] = 0
+        p = lag.Laguerre.fit(x, yw, 3, w=w)
+        assert_almost_equal(p(x), y)
+
+    def test_identity(self) :
+        x = np.linspace(0,3)
+        p = lag.Laguerre.identity()
+        assert_almost_equal(p(x), x)
+        p = lag.Laguerre.identity([1,3])
+        assert_almost_equal(p(x), x)
+#
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/numpy/polynomial/tests/test_legendre.py b/numpy/polynomial/tests/test_legendre.py
index f963bd2df..a35e6605e 100644
--- a/numpy/polynomial/tests/test_legendre.py
+++ b/numpy/polynomial/tests/test_legendre.py
@@ -429,7 +429,7 @@ class TestLegendreClass(TestCase) :
     def test_degree(self) :
         assert_equal(self.p1.degree(), 2)
 
-    def test_trimdeg(self) :
+    def test_cutdeg(self) :
         assert_raises(ValueError, self.p1.cutdeg, .5)
         assert_raises(ValueError, self.p1.cutdeg, -1)
         assert_equal(len(self.p1.cutdeg(3)), 3)
@@ -531,3 +531,7 @@ class TestLegendreClass(TestCase) :
         assert_almost_equal(p(x), x)
         p = leg.Legendre.identity([1,3])
         assert_almost_equal(p(x), x)
+#
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/numpy/polynomial/tests/test_polynomial.py b/numpy/polynomial/tests/test_polynomial.py
index 5890ac13f..4b93ea118 100644
--- a/numpy/polynomial/tests/test_polynomial.py
+++ b/numpy/polynomial/tests/test_polynomial.py
@@ -400,7 +400,7 @@ class TestPolynomialClass(TestCase) :
     def test_degree(self) :
         assert_equal(self.p1.degree(), 2)
 
-    def test_trimdeg(self) :
+    def test_cutdeg(self) :
         assert_raises(ValueError, self.p1.cutdeg, .5)
         assert_raises(ValueError, self.p1.cutdeg, -1)
         assert_equal(len(self.p1.cutdeg(3)), 3)
@@ -502,3 +502,7 @@ class TestPolynomialClass(TestCase) :
         assert_almost_equal(p(x), x)
         p = poly.Polynomial.identity([1,3])
         assert_almost_equal(p(x), x)
+#
+
+if __name__ == "__main__":
+    run_module_suite()