ENH: First commit of hermite and laguerre polynomials. The documentation and

tests still need fixes.

6 files changed, 5022 insertions, 0 deletions

diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
new file mode 100644
index 000000000..550e316de
--- /dev/null
+++ b/numpy/polynomial/hermite.py
@@ -0,0 +1,1141 @@
+"""
+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 import polynomial as P
+    >>> p = P.Polynomial(np.arange(4))
+    >>> p
+    Polynomial([ 0.,  1.,  2.,  3.], [-1.,  1.])
+    >>> c = P.Hermite(P.poly2herm(p.coef))
+    >>> c
+    Hermite([ 1.  ,  3.25,  1.  ,  0.75], [-1.,  1.])
+
+    """
+    [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
+    --------
+    >>> c = P.Hermite(range(4))
+    >>> c
+    Hermite([ 0.,  1.,  2.,  3.], [-1.,  1.])
+    >>> p = c.convert(kind=P.Polynomial)
+    >>> p
+    Polynomial([-1. , -3.5,  3. ,  7.5], [-1.,  1.])
+    >>> P.herm2poly(range(4))
+    array([-1. , -3.5,  3. ,  7.5])
+
+
+    """
+    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
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.hermline(3,2)
+    array([3, 2])
+    >>> L.hermval(-3, L.hermline(3,2)) # should be -3
+    -3.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
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.hermfromroots((-1,0,1)) # x^3 - x relative to the standard basis
+    array([ 0. , -0.4,  0. ,  0.4])
+    >>> j = complex(0,1)
+    >>> L.hermfromroots((-j,j)) # x^2 + 1 relative to the standard basis
+    array([ 1.33333333+0.j,  0.00000000+0.j,  0.66666667+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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.hermadd(c1,c2)
+    array([ 4.,  4.,  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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.hermsub(c1,c2)
+    array([-2.,  0.,  2.])
+    >>> L.hermsub(c2,c1) # -C.hermsub(c1,c2)
+    array([ 2.,  0., -2.])
+
+    """
+    # 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) = ((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]/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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2)
+    >>> P.hermmul(c1,c2) # multiplication requires "reprojection"
+    array([  4.33333333,  10.4       ,  11.66666667,   3.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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.hermdiv(c1,c2) # quotient "intuitive," remainder not
+    (array([ 3.]), array([-8., -4.]))
+    >>> c2 = (0,1,2,3)
+    >>> L.hermdiv(c2,c1) # neither "intuitive"
+    (array([-0.07407407,  1.66666667]), array([-1.03703704, -2.51851852]))
+
+    """
+    # 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
+    --------
+
+    """
+    # 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 import legendre as L
+    >>> cs = (1,2,3,4)
+    >>> L.hermder(cs)
+    array([  6.,   9.,  20.])
+    >>> L.hermder(cs,3)
+    array([ 60.])
+    >>> L.hermder(cs,scl=-1)
+    array([ -6.,  -9., -20.])
+    >>> L.hermder(cs,2,-1)
+    array([  9.,  60.])
+
+    """
+    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 import legendre as L
+    >>> cs = (1,2,3)
+    >>> L.hermint(cs)
+    array([ 0.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.hermint(cs,3)
+    array([  1.66666667e-02,  -1.78571429e-02,   4.76190476e-02,
+            -1.73472348e-18,   1.90476190e-02,   9.52380952e-03])
+    >>> L.hermint(cs, k=3)
+    array([ 3.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.hermint(cs, lbnd=-2)
+    array([ 7.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.hermint(cs, scl=2)
+    array([ 0.66666667,  0.8       ,  1.33333333,  1.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]*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
+    --------
+
+    """
+    # 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.
+
+    """
+    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
+    --------
+
+    """
+    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
+    --------
+    >>> import numpy.polynomial as P
+    >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots
+    array([-0.60582959+0.j        , -0.07208521-0.63832674j,
+           -0.07208521+0.63832674j])
+    >>> P.hermroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots
+    array([-0.85099543, -0.11407192,  0.51506735])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) <= 1 :
+        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..36e452074
--- /dev/null
+++ b/numpy/polynomial/hermite_e.py
@@ -0,0 +1,1138 @@
+"""
+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
+
+__all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
+        'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermeval',
+        'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots',
+        'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'Hermite_e']
+
+import numpy as np
+import numpy.linalg as la
+import polyutils as pu
+import warnings
+from polytemplate import polytemplate
+
+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 import polynomial as P
+    >>> p = P.Polynomial(np.arange(4))
+    >>> p
+    Polynomial([ 0.,  1.,  2.,  3.], [-1.,  1.])
+    >>> c = P.Hermite(P.poly2herme(p.coef))
+    >>> c
+    Hermite([ 1.  ,  3.25,  1.  ,  0.75], [-1.,  1.])
+
+    """
+    [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
+    --------
+    >>> c = P.Hermite(range(4))
+    >>> c
+    Hermite([ 0.,  1.,  2.,  3.], [-1.,  1.])
+    >>> p = c.convert(kind=P.Polynomial)
+    >>> p
+    Polynomial([-1. , -3.5,  3. ,  7.5], [-1.,  1.])
+    >>> P.herme2poly(range(4))
+    array([-1. , -3.5,  3. ,  7.5])
+
+
+    """
+    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
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.hermeline(3,2)
+    array([3, 2])
+    >>> L.hermeval(-3, L.hermeline(3,2)) # should be -3
+    -3.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
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.hermefromroots((-1,0,1)) # x^3 - x relative to the standard basis
+    array([ 0. , -0.4,  0. ,  0.4])
+    >>> j = complex(0,1)
+    >>> L.hermefromroots((-j,j)) # x^2 + 1 relative to the standard basis
+    array([ 1.33333333+0.j,  0.00000000+0.j,  0.66666667+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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.hermeadd(c1,c2)
+    array([ 4.,  4.,  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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.hermesub(c1,c2)
+    array([-2.,  0.,  2.])
+    >>> L.hermesub(c2,c1) # -C.hermesub(c1,c2)
+    array([ 2.,  0., -2.])
+
+    """
+    # 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) = ((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]/2
+    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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2)
+    >>> P.hermemul(c1,c2) # multiplication requires "reprojection"
+    array([  4.33333333,  10.4       ,  11.66666667,   3.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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.hermediv(c1,c2) # quotient "intuitive," remainder not
+    (array([ 3.]), array([-8., -4.]))
+    >>> c2 = (0,1,2,3)
+    >>> L.hermediv(c2,c1) # neither "intuitive"
+    (array([-0.07407407,  1.66666667]), array([-1.03703704, -2.51851852]))
+
+    """
+    # 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
+    --------
+
+    """
+    # 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 import legendre as L
+    >>> cs = (1,2,3,4)
+    >>> L.hermeder(cs)
+    array([  6.,   9.,  20.])
+    >>> L.hermeder(cs,3)
+    array([ 60.])
+    >>> L.hermeder(cs,scl=-1)
+    array([ -6.,  -9., -20.])
+    >>> L.hermeder(cs,2,-1)
+    array([  9.,  60.])
+
+    """
+    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 import legendre as L
+    >>> cs = (1,2,3)
+    >>> L.hermeint(cs)
+    array([ 0.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.hermeint(cs,3)
+    array([  1.66666667e-02,  -1.78571429e-02,   4.76190476e-02,
+            -1.73472348e-18,   1.90476190e-02,   9.52380952e-03])
+    >>> L.hermeint(cs, k=3)
+    array([ 3.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.hermeint(cs, lbnd=-2)
+    array([ 7.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.hermeint(cs, scl=2)
+    array([ 0.66666667,  0.8       ,  1.33333333,  1.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]*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
+    --------
+
+    """
+    # 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.
+
+    """
+    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
+    --------
+
+    """
+    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
+    --------
+    >>> import numpy.polynomial as P
+    >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots
+    array([-0.60582959+0.j        , -0.07208521-0.63832674j,
+           -0.07208521+0.63832674j])
+    >>> P.hermeroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots
+    array([-0.85099543, -0.11407192,  0.51506735])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) <= 1 :
+        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
+
+
+#
+# Hermite_e series class
+#
+
+exec polytemplate.substitute(name='Hermite_e', nick='herme', domain='[-1,1]')
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
new file mode 100644
index 000000000..94c495deb
--- /dev/null
+++ b/numpy/polynomial/laguerre.py
@@ -0,0 +1,1141 @@
+"""
+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 import polynomial as P
+    >>> p = P.Polynomial(np.arange(4))
+    >>> p
+    Polynomial([ 0.,  1.,  2.,  3.], [-1.,  1.])
+    >>> c = P.Laguerre(P.poly2lag(p.coef))
+    >>> c
+    Laguerre([ 1.  ,  3.25,  1.  ,  0.75], [-1.,  1.])
+
+    """
+    [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
+    --------
+    >>> c = P.Laguerre(range(4))
+    >>> c
+    Laguerre([ 0.,  1.,  2.,  3.], [-1.,  1.])
+    >>> p = c.convert(kind=P.Polynomial)
+    >>> p
+    Polynomial([-1. , -3.5,  3. ,  7.5], [-1.,  1.])
+    >>> P.lag2poly(range(4))
+    array([-1. , -3.5,  3. ,  7.5])
+
+
+    """
+    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
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.lagline(3,2)
+    array([3, 2])
+    >>> L.lagval(-3, L.lagline(3,2)) # should be -3
+    -3.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
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.lagfromroots((-1,0,1)) # x^3 - x relative to the standard basis
+    array([ 0. , -0.4,  0. ,  0.4])
+    >>> j = complex(0,1)
+    >>> L.lagfromroots((-j,j)) # x^2 + 1 relative to the standard basis
+    array([ 1.33333333+0.j,  0.00000000+0.j,  0.66666667+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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.lagadd(c1,c2)
+    array([ 4.,  4.,  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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.lagsub(c1,c2)
+    array([-2.,  0.,  2.])
+    >>> L.lagsub(c2,c1) # -C.lagsub(c1,c2)
+    array([ 2.,  0., -2.])
+
+    """
+    # 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) + 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]
+    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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2)
+    >>> P.lagmul(c1,c2) # multiplication requires "reprojection"
+    array([  4.33333333,  10.4       ,  11.66666667,   3.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 = 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 import legendre as L
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2,1)
+    >>> L.lagdiv(c1,c2) # quotient "intuitive," remainder not
+    (array([ 3.]), array([-8., -4.]))
+    >>> c2 = (0,1,2,3)
+    >>> L.lagdiv(c2,c1) # neither "intuitive"
+    (array([-0.07407407,  1.66666667]), array([-1.03703704, -2.51851852]))
+
+    """
+    # 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
+    --------
+
+    """
+    # 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 import legendre as L
+    >>> cs = (1,2,3,4)
+    >>> L.lagder(cs)
+    array([  6.,   9.,  20.])
+    >>> L.lagder(cs,3)
+    array([ 60.])
+    >>> L.lagder(cs,scl=-1)
+    array([ -6.,  -9., -20.])
+    >>> L.lagder(cs,2,-1)
+    array([  9.,  60.])
+
+    """
+    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 import legendre as L
+    >>> cs = (1,2,3)
+    >>> L.lagint(cs)
+    array([ 0.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.lagint(cs,3)
+    array([  1.66666667e-02,  -1.78571429e-02,   4.76190476e-02,
+            -1.73472348e-18,   1.90476190e-02,   9.52380952e-03])
+    >>> L.lagint(cs, k=3)
+    array([ 3.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.lagint(cs, lbnd=-2)
+    array([ 7.33333333,  0.4       ,  0.66666667,  0.6       ])
+    >>> L.lagint(cs, scl=2)
+    array([ 0.66666667,  0.8       ,  1.33333333,  1.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
+    --------
+
+    """
+    # 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.
+
+    """
+    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
+    --------
+
+    """
+    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
+    --------
+    >>> import numpy.polynomial as P
+    >>> P.polyroots((1, 2, 3, 4)) # 4x^3 + 3x^2 + 2x + 1 has two complex roots
+    array([-0.60582959+0.j        , -0.07208521-0.63832674j,
+           -0.07208521+0.63832674j])
+    >>> P.lagroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0 has only real roots
+    array([-0.85099543, -0.11407192,  0.51506735])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) <= 1 :
+        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/tests/test_hermite.py b/numpy/polynomial/tests/test_hermite.py
new file mode 100644
index 000000000..7e361e804
--- /dev/null
+++ b/numpy/polynomial/tests/test_hermite.py
@@ -0,0 +1,536 @@
+"""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_trimdeg(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..aa01baf8e
--- /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,.5])
+        for i in range(1, 5):
+            ser = [0]*i + [1]
+            tgt = [0]*(i - 1) + [i, 0, .5]
+            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, .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]
+            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]), Hlist[i])
+
+    def test_poly2herme(self) :
+        for i in range(10) :
+            assert_almost_equal(herme.poly2herme(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 TestHermite_eClass(TestCase) :
+
+    p1 = herme.Hermite_e([1,2,3])
+    p2 = herme.Hermite_e([1,2,3], [0,1])
+    p3 = herme.Hermite_e([1,2])
+    p4 = herme.Hermite_e([2,2,3])
+    p5 = herme.Hermite_e([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.Hermite_e([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.Hermite_e([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.Hermite_e([ 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 = herme.Hermite_e([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.Hermite_e([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.Hermite_e([1])
+        trem = herme.Hermite_e([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.Hermite_e([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_trimdeg(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.Hermite_e(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.Hermite_e(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.Hermite_e.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.Hermite_e.fit(x, y, 3)
+        assert_almost_equal(p.domain, [0,3])
+
+        # test that fit works in given domains
+        p = herme.Hermite_e.fit(x, y, 3, None)
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [0,3])
+        p = herme.Hermite_e.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.Hermite_e.fit(x, yw, 3, w=w)
+        assert_almost_equal(p(x), y)
+
+    def test_identity(self) :
+        x = np.linspace(0,3)
+        p = herme.Hermite_e.identity()
+        assert_almost_equal(p(x), x)
+        p = herme.Hermite_e.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..4d1b80e4b
--- /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_trimdeg(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()