ENH: Add support for Legendre polynomials.

5 files changed, 1788 insertions, 4 deletions

diff --git a/numpy/polynomial/__init__.py b/numpy/polynomial/__init__.py
index 7e755ca52..6b004acc4 100644
--- a/numpy/polynomial/__init__.py
+++ b/numpy/polynomial/__init__.py
@@ -15,6 +15,7 @@ information can be found in the docstring for the module of interest.
 """
 from polynomial import *
 from chebyshev import *
+from legendre import *
 from polyutils import *
 
 from numpy.testing import Tester
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index 23f8a728f..ea064a695 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -408,9 +408,10 @@ def cheb2poly(cs) :
     else:
         c0 = cs[-2]
         c1 = cs[-1]
-        for i in range(n - 3, -1, -1) :
+        # i is the current degree of c1
+        for i in range(n - 1, 1, -1) :
             tmp = c0
-            c0 = polysub(cs[i], c1)
+            c0 = polysub(cs[i - 2], c1)
             c1 = polyadd(tmp, polymulx(c1)*2)
         return polyadd(c0, polymulx(c1))
 
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
new file mode 100644
index 000000000..acd2d0554
--- /dev/null
+++ b/numpy/polynomial/legendre.py
@@ -0,0 +1,1249 @@
+"""
+Objects for dealing with Legendre series.
+
+This module provides a number of objects (mostly functions) useful for
+dealing with Legendre series, including a `Legendre` class that
+encapsulates the usual arithmetic operations.  (General information
+on how this module represents and works with such polynomials is in the
+docstring for its "parent" sub-package, `numpy.polynomial`).
+
+Constants
+---------
+- `legdomain` -- Legendre series default domain, [-1,1].
+- `legzero` -- Legendre series that evaluates identically to 0.
+- `legone` -- Legendre series that evaluates identically to 1.
+- `legx` -- Legendre series for the identity map, ``f(x) = x``.
+
+Arithmetic
+----------
+- `legmulx` -- multiply a Legendre series in ``P_i(x)`` by ``x``.
+- `legadd` -- add two Legendre series.
+- `legsub` -- subtract one Legendre series from another.
+- `legmul` -- multiply two Legendre series.
+- `legdiv` -- divide one Legendre series by another.
+- `legval` -- evaluate a Legendre series at given points.
+
+Calculus
+--------
+- `legder` -- differentiate a Legendre series.
+- `legint` -- integrate a Legendre series.
+
+Misc Functions
+--------------
+- `legfromroots` -- create a Legendre series with specified roots.
+- `legroots` -- find the roots of a Legendre series.
+- `legvander` -- Vandermonde-like matrix for Legendre polynomials.
+- `legfit` -- least-squares fit returning a Legendre series.
+- `legtrim` -- trim leading coefficients from a Legendre series.
+- `legline` -- Legendre series of given straight line.
+- `leg2poly` -- convert a Legendre series to a polynomial.
+- `poly2leg` -- convert a polynomial to a Legendre series.
+
+Classes
+-------
+- `Legendre` -- A Legendre series class.
+
+See also
+--------
+`numpy.polynomial`
+
+Notes
+-----
+The implementations of multiplication, division, integration, and
+differentiation use the algebraic identities [1]_:
+
+.. math ::
+    T_n(x) = \\frac{z^n + z^{-n}}{2} \\\\
+    z\\frac{dx}{dz} = \\frac{z - z^{-1}}{2}.
+
+where
+
+.. math :: x = \\frac{z + z^{-1}}{2}.
+
+These identities allow a Chebyshev series to be expressed as a finite,
+symmetric Laurent series.  In this module, this sort of Laurent series
+is referred to as a "z-series."
+
+References
+----------
+.. [1] A. T. Benjamin, et al., "Combinatorial Trigonometry with Chebyshev
+  Polynomials," *Journal of Statistical Planning and Inference 14*, 2008
+  (preprint: http://www.math.hmc.edu/~benjamin/papers/CombTrig.pdf, pg. 4)
+
+"""
+from __future__ import division
+
+__all__ = ['legzero', 'legone', 'legx', 'legdomain', 'legline',
+        'legadd', 'legsub', 'legmulx', 'legmul', 'legdiv', 'legval',
+        'legder', 'legint', 'leg2poly', 'poly2leg', 'legfromroots',
+        'legvander', 'legfit', 'legtrim', 'legroots', 'Legendre']
+
+import numpy as np
+import numpy.linalg as la
+import polyutils as pu
+import warnings
+from polytemplate import polytemplate
+
+legtrim = pu.trimcoef
+
+def poly2leg(pol) :
+    """
+    poly2leg(pol)
+
+    Convert a polynomial to a Legendre 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 Legendre 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 Legendre
+        series.
+
+    See Also
+    --------
+    leg2poly
+
+    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.Legendre(P.poly2leg(p.coef))
+    >>> c
+    Legendre([ 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 = legadd(legmulx(res), pol[i])
+    return res
+
+
+def leg2poly(cs) :
+    """
+    Convert a Legendre series to a polynomial.
+
+    Convert an array representing the coefficients of a Legendre 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 Legendre 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
+    --------
+    poly2leg
+
+    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
+    >>> c = P.Chebyshev(np.arange(4))
+    >>> c
+    Chebyshev([ 0.,  1.,  2.,  3.], [-1.,  1.])
+    >>> p = P.Polynomial(P.cheb2poly(c.coef))
+    >>> p
+    Polynomial([ -2.,  -8.,   4.,  12.], [-1.,  1.])
+
+    """
+    from polynomial import polyadd, polysub, polymulx
+
+    [cs] = pu.as_series([cs])
+    n = len(cs)
+    if n < 3:
+        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, (polymulx(c1)*(2*i - 1))/i)
+        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.
+#
+
+# Legendre
+legdomain = np.array([-1,1])
+
+# Legendre coefficients representing zero.
+legzero = np.array([0])
+
+# Legendre coefficients representing one.
+legone = np.array([1])
+
+# Legendre coefficients representing the identity x.
+legx = np.array([0,1])
+
+
+def legline(off, scl) :
+    """
+    Legendre 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 Legendre series for
+        ``off + scl*x``.
+
+    See Also
+    --------
+    polyline, chebline
+
+    Examples
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.legline(3,2)
+    array([3, 2])
+    >>> L.legval(-3, L.chebline(3,2)) # should be -3
+    -3.0
+
+    """
+    if scl != 0 :
+        return np.array([off,scl])
+    else :
+        return np.array([off])
+
+
+def legtimesx(cs):
+    """Multiply a Legendre series by x.
+
+    Multiply the Legendre series `cs` by x, where x is the independent
+    variable.
+
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Legendre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the result of the multiplication.
+
+    Notes
+    -----
+    The multiplication uses the recursion relationship for Legendre
+    polynomials in the form
+
+    .. math::
+
+    xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1)
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    # The zero series needs special treatment
+    if len(cs) == 1 and cs[0] == 0:
+        return cs
+
+    prd = np.empty(len(cs) + 1, dtype=cs.dtype)
+    prd[0] = cs[0]*0
+    prd[1] = cs[0]
+    for i in range(1, len(cs)):
+        j = i + 1
+        k = i - 1
+        s = i + j
+        prd[j] = (cs[i]*j)/s
+        prd[k] += (cs[i]*i)/s
+    return prd
+
+
+def chebline(off, scl) :
+    """
+    Chebyshev 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 Chebyshev series for
+        ``off + scl*x``.
+
+    See Also
+    --------
+    polyline
+
+    Examples
+    --------
+    >>> import numpy.polynomial.chebyshev as C
+    >>> C.chebline(3,2)
+    array([3, 2])
+    >>> C.chebval(-3, C.chebline(3,2)) # should be -3
+    -3.0
+
+    """
+    if scl != 0 :
+        return np.array([off,scl])
+    else :
+        return np.array([off])
+
+
+def legfromroots(roots) :
+    """
+    Generate a Legendre 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 Legendre 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 Legendre
+    (basis) polynomial over the domain `[-1,1]`.  Note that, unlike
+    `polyfromroots`, due to the nature of the Legendre basis set, the
+    above identity *does not* imply :math:`c_n = 1` identically (see
+    Examples).
+
+    Examples
+    --------
+    >>> import numpy.polynomial.legendre as L
+    >>> L.legfromroots((-1,0,1)) # x^3 - x relative to the standard basis
+    array([ 0. , -0.4,  0. ,  0.4])
+    >>> j = complex(0,1)
+    >>> L.legfromroots((-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 = legsub(legmulx(prd), r*prd)
+        return prd
+
+
+def legadd(c1, c2):
+    """
+    Add one Legendre series to another.
+
+    Returns the sum of two Legendre 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 Legendre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the Legendre series of their sum.
+
+    See Also
+    --------
+    legsub, legmul, legdiv, legpow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the sum of two Legendre series
+    is a Legendre 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.legadd(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 legsub(c1, c2):
+    """
+    Subtract one Legendre series from another.
+
+    Returns the difference of two Legendre 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 Legendre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Legendre series coefficients representing their difference.
+
+    See Also
+    --------
+    legadd, legmul, legdiv, legpow
+
+    Notes
+    -----
+    Unlike multiplication, division, etc., the difference of two Legendre
+    series is a Legendre 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.legsub(c1,c2)
+    array([-2.,  0.,  2.])
+    >>> L.legsub(c2,c1) # -C.legsub(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 legmulx(cs):
+    """Multiply a Legendre series by x.
+
+    Multiply the Legendre series `cs` by x, where x is the independent
+    variable.
+
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Legendre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Array representing the result of the multiplication.
+
+    Notes
+    -----
+    The multiplication uses the recursion relationship for Legendre
+    polynomials in the form
+
+    .. math::
+
+    xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1)
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    # The zero series needs special treatment
+    if len(cs) == 1 and cs[0] == 0:
+        return cs
+
+    prd = np.empty(len(cs) + 1, dtype=cs.dtype)
+    prd[0] = cs[0]*0
+    prd[1] = cs[0]
+    for i in range(1, len(cs)):
+        j = i + 1
+        k = i - 1
+        s = i + j
+        prd[j] = (cs[i]*j)/s
+        prd[k] += (cs[i]*i)/s
+    return prd
+
+
+def legmul(c1, c2):
+    """
+    Multiply one Legendre series by another.
+
+    Returns the product of two Legendre 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 Legendre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    out : ndarray
+        Of Legendre series coefficients representing their product.
+
+    See Also
+    --------
+    legadd, legsub, legdiv, legpow
+
+    Notes
+    -----
+    In general, the (polynomial) product of two C-series results in terms
+    that are not in the Chebyshev polynomial basis set.  Thus, to express
+    the product as a C-series, it is typically necessary to "re-project"
+    the product onto said basis set, which typically produces
+    "un-intuitive" (but correct) results; see Examples section below.
+
+    Examples
+    --------
+    >>> from numpy.polynomial import legendre as P
+    >>> c1 = (1,2,3)
+    >>> c2 = (3,2)
+    >>> P.legmul(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 = legsub(cs[-i]*xs, (c1*(nd - 1))/nd)
+            c1 = legadd(tmp, (legmulx(c1)*(2*nd - 1))/nd)
+    return legadd(c0, legmulx(c1))
+
+
+def legdiv(c1, c2):
+    """
+    Divide one Legendre series by another.
+
+    Returns the quotient-with-remainder of two Legendre 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 Legendre series coefficients ordered from low to
+        high.
+
+    Returns
+    -------
+    [quo, rem] : ndarrays
+        Of Legendre series coefficients representing the quotient and
+        remainder.
+
+    See Also
+    --------
+    legadd, legsub, legmul, legpow
+
+    Notes
+    -----
+    In general, the (polynomial) division of one Legendre series by another
+    results in quotient and remainder terms that are not in the Legendre
+    polynomial basis set.  Thus, to express these results as a Legendre
+    series, it is necessary to "re-project" the results onto the Legendre
+    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.legdiv(c1,c2) # quotient "intuitive," remainder not
+    (array([ 3.]), array([-8., -4.]))
+    >>> c2 = (0,1,2,3)
+    >>> L.legdiv(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 = legmul([0]*i + [1], c2)
+            q = rem[-1]/p[-1]
+            rem = rem[:-1] - q*p[:-1]
+            quo[i] = q
+        return quo, pu.trimseq(rem)
+
+
+def legpow(cs, pow, maxpower=16) :
+    """Raise a Legendre series to a power.
+
+    Returns the Legendre 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 Legendre 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
+        Legendre series of power.
+
+    See Also
+    --------
+    legadd, legsub, legmul, legdiv
+
+    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 = legmul(prd, cs)
+        return prd
+
+
+def legder(cs, m=1, scl=1) :
+    """
+    Differentiate a Legendre 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 Legendre 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
+        Legendre series of the derivative.
+
+    See Also
+    --------
+    legint
+
+    Notes
+    -----
+    In general, the result of differentiating a Legendre 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.legder(cs)
+    array([  6.,   9.,  20.])
+    >>> L.legder(cs,3)
+    array([ 60.])
+    >>> L.legder(cs,scl=-1)
+    array([ -6.,  -9., -20.])
+    >>> L.legder(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 - 1)*cs[j]
+                cs[j - 2] += cs[j]
+            cs = der
+        return cs
+
+
+def legint(cs, m=1, k=[], lbnd=0, scl=1):
+    """
+    Integrate a Legendre series.
+
+    Returns a Legendre series that is the Legendre 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 Legendre 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 Legendre 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
+        Legendre series coefficients of the integral.
+
+    Raises
+    ------
+    ValueError
+        If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or
+        ``np.isscalar(scl) == False``.
+
+    See Also
+    --------
+    legder
+
+    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 legyshev as L
+    >>> cs = (1,2,3)
+    >>> L.legint(cs)
+    array([ 0.5, -0.5,  0.5,  0.5])
+    >>> L.legint(cs,3)
+    array([ 0.03125   , -0.1875    ,  0.04166667, -0.05208333,  0.01041667,
+            0.00625   ])
+    >>> L.legint(cs, k=3)
+    array([ 3.5, -0.5,  0.5,  0.5])
+    >>> L.legint(cs,lbnd=-2)
+    array([ 8.5, -0.5,  0.5,  0.5])
+    >>> L.legint(cs,scl=-2)
+    array([-1.,  1., -1., -1.])
+
+    """
+    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):
+                t = cs[j]/(2*j + 1)
+                tmp[j + 1] = t
+                tmp[j - 1] -= t
+            tmp[0] += k[i] - legval(lbnd, tmp)
+            cs = tmp
+    return cs
+
+
+def legval(x, cs):
+    """Evaluate a Legendre 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 Chebyshev 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
+    --------
+    legfit
+
+    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*x*(2*nd - 1))/nd
+    return c0 + c1*x
+
+
+def legvander(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 Legendre 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 Legendre 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)
+    # Use forward recursion to generate the entries. This is not as accurate
+    # as reverse recursion in this application but it is more efficient.
+    v[0] = x*0 + 1
+    if ideg > 0 :
+        v[1] = x
+        for i in range(2, ideg + 1) :
+            v[i] = (v[i-1]*x*(2*i - 1) - v[i-2]*(i - 1))/i
+    return np.rollaxis(v, 0, v.ndim)
+
+
+def legfit(x, y, deg, rcond=None, full=False, w=None):
+    """
+    Least squares fit of Legendre series to data.
+
+    Fit a Legendre series ``p(x) = p[0] * P_{0}(x) + ... + p[deg] *
+    P_{deg}(x)`` of degree `deg` to points `(x, y)`. Returns a vector of
+    coefficients `p` that minimises the squared error.
+
+    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)
+        Legendre 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
+    --------
+    legval : Evaluates a Legendre series.
+    legvander : Vandermonde matrix of Legendre 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 Legendre 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 Legendre 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 = legvander(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 legroots(cs):
+    """
+    Compute the roots of a Chebyshev series.
+
+    Return the roots (a.k.a "zeros") of the Legendre series represented by
+    `cs`, which is the sequence of the C-series' coefficients from lowest
+    order "term" to highest, e.g., [1,2,3] represents the Legendre series
+    ``P_0 + 2*P_1 + 3*P_2``.
+
+    Parameters
+    ----------
+    cs : array_like
+        1-d array of Legendre 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 Legendre 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
+    >>> import numpy.polynomial.Legendre as L
+    >>> P.polyroots((-1,1,-1,1)) # x^3 - x^2 + x - 1 has two complex roots
+    array([ -4.99600361e-16-1.j,  -4.99600361e-16+1.j,   1.00000e+00+0.j])
+    >>> L.legroots((-1,1,-1,1)) # T3 - T2 + T1 - T0 has only real roots
+    array([ -5.00000000e-01,   2.60860684e-17,   1.00000000e+00])
+
+    """
+    # cs is a trimmed copy
+    [cs] = pu.as_series([cs])
+    if len(cs) <= 1 :
+        return np.array([], dtype=cs.dtype)
+    if len(cs) == 2 :
+        return np.array([-cs[0]/cs[1]])
+
+    n = len(cs) - 1
+    cs /= cs[-1]
+    cmat = np.zeros((n,n), dtype=cs.dtype)
+    cmat[1, 0] = 1
+    for i in range(1, n):
+        tmp = 2*i + 1
+        cmat[i - 1, i] = i/tmp
+        if i != n - 1:
+            cmat[i + 1, i] = (i + 1)/tmp
+        else:
+            cmat[:, i] -= cs[:-1]*(i + 1)/tmp
+    roots = la.eigvals(cmat)
+    roots.sort()
+    return roots
+
+
+#
+# Legendre series class
+#
+
+exec polytemplate.substitute(name='Legendre', nick='leg', domain='[-1,1]')
diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py
index ffa8e984d..65bb877f4 100644
--- a/numpy/polynomial/tests/test_chebyshev.py
+++ b/numpy/polynomial/tests/test_chebyshev.py
@@ -339,11 +339,11 @@ class TestMisc(TestCase) :
 
     def test_cheb2poly(self) :
         for i in range(10) :
-            assert_equal(ch.cheb2poly([0]*i + [1]), Tlist[i])
+            assert_almost_equal(ch.cheb2poly([0]*i + [1]), Tlist[i])
 
     def test_poly2cheb(self) :
         for i in range(10) :
-            assert_equal(ch.poly2cheb(Tlist[i]), [0]*i + [1])
+            assert_almost_equal(ch.poly2cheb(Tlist[i]), [0]*i + [1])
 
     def test_chebpts1(self):
         #test exceptions
diff --git a/numpy/polynomial/tests/test_legendre.py b/numpy/polynomial/tests/test_legendre.py
new file mode 100644
index 000000000..f963bd2df
--- /dev/null
+++ b/numpy/polynomial/tests/test_legendre.py
@@ -0,0 +1,533 @@
+"""Tests for legendre module.
+
+"""
+from __future__ import division
+
+import numpy as np
+import numpy.polynomial.legendre as leg
+import numpy.polynomial.polynomial as poly
+from numpy.testing import *
+
+P0 = np.array([ 1])
+P1 = np.array([ 0,  1])
+P2 = np.array([-1,  0,    3])/2
+P3 = np.array([ 0, -3,    0,    5])/2
+P4 = np.array([ 3,  0,  -30,    0,  35])/8
+P5 = np.array([ 0, 15,    0,  -70,   0,   63])/8
+P6 = np.array([-5,  0,  105,    0,-315,    0,   231])/16
+P7 = np.array([ 0,-35,    0,  315,   0, -693,     0,   429])/16
+P8 = np.array([35,  0,-1260,    0,6930,    0,-12012,     0,6435])/128
+P9 = np.array([ 0,315,    0,-4620,   0,18018,     0,-25740,   0,12155])/128
+
+Plist = [P0, P1, P2, P3, P4, P5, P6, P7, P8, P9]
+
+def trim(x) :
+    return leg.legtrim(x, tol=1e-6)
+
+
+class TestConstants(TestCase) :
+
+    def test_legdomain(self) :
+        assert_equal(leg.legdomain, [-1, 1])
+
+    def test_legzero(self) :
+        assert_equal(leg.legzero, [0])
+
+    def test_legone(self) :
+        assert_equal(leg.legone, [1])
+
+    def test_legx(self) :
+        assert_equal(leg.legx, [0, 1])
+
+
+class TestArithmetic(TestCase) :
+    x = np.linspace(-1, 1, 100)
+    y0 = poly.polyval(x, P0)
+    y1 = poly.polyval(x, P1)
+    y2 = poly.polyval(x, P2)
+    y3 = poly.polyval(x, P3)
+    y4 = poly.polyval(x, P4)
+    y5 = poly.polyval(x, P5)
+    y6 = poly.polyval(x, P6)
+    y7 = poly.polyval(x, P7)
+    y8 = poly.polyval(x, P8)
+    y9 = poly.polyval(x, P9)
+    y = [y0, y1, y2, y3, y4, y5, y6, y7, y8, y9]
+
+    def test_legval(self) :
+        def f(x) :
+            return x*(x**2 - 1)
+
+        #check empty input
+        assert_equal(leg.legval([], [1]).size, 0)
+
+        #check normal input)
+        for i in range(10) :
+            msg = "At i=%d" % i
+            ser = np.zeros
+            tgt = self.y[i]
+            res = leg.legval(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(leg.legval(x, [1]).shape, dims)
+            assert_equal(leg.legval(x, [1,0]).shape, dims)
+            assert_equal(leg.legval(x, [1,0,0]).shape, dims)
+
+    def test_legadd(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 = leg.legadd([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_legsub(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 = leg.legsub([0]*i + [1], [0]*j + [1])
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+    def test_legmulx(self):
+        assert_equal(leg.legmulx([0]), [0])
+        assert_equal(leg.legmulx([1]), [0,1])
+        for i in range(1, 5):
+            tmp = 2*i + 1
+            ser = [0]*i + [1]
+            tgt = [0]*(i - 1) + [i/tmp, 0, (i + 1)/tmp]
+            assert_equal(leg.legmulx(ser), tgt)
+
+    def test_legmul(self) :
+        # check values of result
+        for i in range(5) :
+            pol1 = [0]*i + [1]
+            val1 = leg.legval(self.x, pol1)
+            for j in range(5) :
+                msg = "At i=%d, j=%d" % (i,j)
+                pol2 = [0]*j + [1]
+                val2 = leg.legval(self.x, pol2)
+                pol3 = leg.legmul(pol1, pol2)
+                val3 = leg.legval(self.x, pol3)
+                assert_(len(pol3) == i + j + 1, msg)
+                assert_almost_equal(val3, val1*val2, err_msg=msg)
+
+    def test_legdiv(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 = leg.legadd(ci, cj)
+                quo, rem = leg.legdiv(tgt, ci)
+                res = leg.legadd(leg.legmul(quo, ci), rem)
+                assert_equal(trim(res), trim(tgt), err_msg=msg)
+
+
+class TestCalculus(TestCase) :
+
+    def test_legint(self) :
+        # check exceptions
+        assert_raises(ValueError, leg.legint, [0], .5)
+        assert_raises(ValueError, leg.legint, [0], -1)
+        assert_raises(ValueError, leg.legint, [0], 1, [0,0])
+
+        # test integration of zero polynomial
+        for i in range(2, 5):
+            k = [0]*(i - 2) + [1]
+            res = leg.legint([0], m=i, k=k)
+            assert_almost_equal(res, [0, 1])
+
+        # check single integration with integration constant
+        for i in range(5) :
+            scl = i + 1
+            pol = [0]*i + [1]
+            tgt = [i] + [0]*i + [1/scl]
+            legpol = leg.poly2leg(pol)
+            legint = leg.legint(legpol, m=1, k=[i])
+            res = leg.leg2poly(legint)
+            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]
+            legpol = leg.poly2leg(pol)
+            legint = leg.legint(legpol, m=1, k=[i], lbnd=-1)
+            assert_almost_equal(leg.legval(-1, legint), 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]
+            legpol = leg.poly2leg(pol)
+            legint = leg.legint(legpol, m=1, k=[i], scl=2)
+            res = leg.leg2poly(legint)
+            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 = leg.legint(tgt, m=1)
+                res = leg.legint(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 = leg.legint(tgt, m=1, k=[k])
+                res = leg.legint(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 = leg.legint(tgt, m=1, k=[k], lbnd=-1)
+                res = leg.legint(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 = leg.legint(tgt, m=1, k=[k], scl=2)
+                res = leg.legint(pol, m=j, k=range(j), scl=2)
+                assert_almost_equal(trim(res), trim(tgt))
+
+    def test_legder(self) :
+        # check exceptions
+        assert_raises(ValueError, leg.legder, [0], .5)
+        assert_raises(ValueError, leg.legder, [0], -1)
+
+        # check that zeroth deriviative does nothing
+        for i in range(5) :
+            tgt = [1] + [0]*i
+            res = leg.legder(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 = leg.legder(leg.legint(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 = leg.legder(leg.legint(tgt, m=j, scl=2), m=j, scl=.5)
+                assert_almost_equal(trim(res), trim(tgt))
+
+
+class TestMisc(TestCase) :
+
+    def test_legfromroots(self) :
+        res = leg.legfromroots([])
+        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 = leg.legfromroots(roots)
+            res = leg.legval(roots, pol)
+            tgt = 0
+            assert_(len(pol) == i + 1)
+            assert_almost_equal(leg.leg2poly(pol)[-1], 1)
+            assert_almost_equal(res, tgt)
+
+    def test_legroots(self) :
+        assert_almost_equal(leg.legroots([1]), [])
+        assert_almost_equal(leg.legroots([1, 2]), [-.5])
+        for i in range(2,5) :
+            tgt = np.linspace(-1, 1, i)
+            res = leg.legroots(leg.legfromroots(tgt))
+            assert_almost_equal(trim(res), trim(tgt))
+
+    def test_legvander(self) :
+        # check for 1d x
+        x = np.arange(3)
+        v = leg.legvander(x, 3)
+        assert_(v.shape == (3,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], leg.legval(x, coef))
+
+        # check for 2d x
+        x = np.array([[1,2],[3,4],[5,6]])
+        v = leg.legvander(x, 3)
+        assert_(v.shape == (3,2,4))
+        for i in range(4) :
+            coef = [0]*i + [1]
+            assert_almost_equal(v[...,i], leg.legval(x, coef))
+
+    def test_legfit(self) :
+        def f(x) :
+            return x*(x - 1)*(x - 2)
+
+        # Test exceptions
+        assert_raises(ValueError, leg.legfit, [1],    [1],     -1)
+        assert_raises(TypeError,  leg.legfit, [[1]],  [1],      0)
+        assert_raises(TypeError,  leg.legfit, [],     [1],      0)
+        assert_raises(TypeError,  leg.legfit, [1],    [[[1]]],  0)
+        assert_raises(TypeError,  leg.legfit, [1, 2], [1],      0)
+        assert_raises(TypeError,  leg.legfit, [1],    [1, 2],   0)
+        assert_raises(TypeError,  leg.legfit, [1],    [1],   0, w=[[1]])
+        assert_raises(TypeError,  leg.legfit, [1],    [1],   0, w=[1,1])
+
+        # Test fit
+        x = np.linspace(0,2)
+        y = f(x)
+        #
+        coef3 = leg.legfit(x, y, 3)
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(leg.legval(x, coef3), y)
+        #
+        coef4 = leg.legfit(x, y, 4)
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(leg.legval(x, coef4), y)
+        #
+        coef2d = leg.legfit(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 = leg.legfit(x, yw, 3, w=w)
+        assert_almost_equal(wcoef3, coef3)
+        #
+        wcoef2d = leg.legfit(x, np.array([yw,yw]).T, 3, w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3,coef3]).T)
+
+    def test_legtrim(self) :
+        coef = [2, -1, 1, 0]
+
+        # Test exceptions
+        assert_raises(ValueError, leg.legtrim, coef, -1)
+
+        # Test results
+        assert_equal(leg.legtrim(coef), coef[:-1])
+        assert_equal(leg.legtrim(coef, 1), coef[:-3])
+        assert_equal(leg.legtrim(coef, 2), [0])
+
+    def test_legline(self) :
+        assert_equal(leg.legline(3,4), [3, 4])
+
+    def test_leg2poly(self) :
+        for i in range(10) :
+            assert_almost_equal(leg.leg2poly([0]*i + [1]), Plist[i])
+
+    def test_poly2leg(self) :
+        for i in range(10) :
+            assert_almost_equal(leg.poly2leg(Plist[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 TestLegendreClass(TestCase) :
+
+    p1 = leg.Legendre([1,2,3])
+    p2 = leg.Legendre([1,2,3], [0,1])
+    p3 = leg.Legendre([1,2])
+    p4 = leg.Legendre([2,2,3])
+    p5 = leg.Legendre([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 = leg.Legendre([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 = leg.Legendre([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 = leg.Legendre([4.13333333, 8.8, 11.23809524, 7.2, 4.62857143])
+        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 = leg.Legendre([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 = leg.Legendre([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 = leg.Legendre([1])
+        trem = leg.Legendre([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 = leg.Legendre([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*(1.5*x**2 - .5) + 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 = leg.Legendre(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, leg.legint([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, leg.legint([1,2,3], 1, 1, scl=.5))
+        p = self.p2.integ(2, [1, 2])
+        assert_almost_equal(p.coef, leg.legint([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 = leg.Legendre(leg.poly2leg([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 = leg.Legendre.fromroots(roots, domain=[0, 1])
+        res = p.coef
+        tgt = leg.poly2leg([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 = leg.Legendre.fit(x, y, 3)
+        assert_almost_equal(p.domain, [0,3])
+
+        # test that fit works in given domains
+        p = leg.Legendre.fit(x, y, 3, None)
+        assert_almost_equal(p(x), y)
+        assert_almost_equal(p.domain, [0,3])
+        p = leg.Legendre.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 = leg.Legendre.fit(x, yw, 3, w=w)
+        assert_almost_equal(p(x), y)
+
+    def test_identity(self) :
+        x = np.linspace(0,3)
+        p = leg.Legendre.identity()
+        assert_almost_equal(p(x), x)
+        p = leg.Legendre.identity([1,3])
+        assert_almost_equal(p(x), x)