From 8e8b85d404e8453e469b3273528c4dc5f30b819f Mon Sep 17 00:00:00 2001
From: Charles Harris <charlesr.harris@gmail.com>
Date: Fri, 16 Dec 2011 20:34:01 -0700
Subject: ENH: Make derivatives and integrals work on multidimensional array.

An axis keyword was added to the function signatures of xxxder and
xxxint, where xxx is any of poly, cheb, leg, lag, herm, herme. The
evaluation method for the Chebeshev series was also changed to avoid
using z_series and to more closely resemble the other implementations.

At some point the z_series will be removed from the chebyshev module
and only used for trigonometric series.
---
 numpy/polynomial/chebyshev.py  | 160 +++++++++++++++++++++++++--------------
 numpy/polynomial/hermite.py    | 129 ++++++++++++++++++++------------
 numpy/polynomial/hermite_e.py  | 135 ++++++++++++++++++++-------------
 numpy/polynomial/laguerre.py   | 139 +++++++++++++++++++++-------------
 numpy/polynomial/legendre.py   | 166 +++++++++++++++++++++++++----------------
 numpy/polynomial/polynomial.py | 160 ++++++++++++++++++++++++---------------
 6 files changed, 556 insertions(+), 333 deletions(-)

(limited to 'numpy/polynomial')

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index d6ccf25ca..c4d0161a2 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -854,26 +854,33 @@ def chebpow(cs, pow, maxpower=16) :
         return _zseries_to_cseries(prd)
 
 
-def chebder(cs, m=1, scl=1) :
+def chebder(c, m=1, scl=1, axis=0) :
     """
     Differentiate a Chebyshev 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
-    ``T_0 + 2*T_1 + 3*T_2``.
+    Returns the Chebyshev series coefficients `c` differentiated `m` times
+    along `axis`.  At each iteration the result is multiplied by `scl` (the
+    scaling factor is for use in a linear change of variable). The argument
+    `c` is an array of coefficients from low to high degree along each
+    axis, e.g., [1,2,3] represents the series ``1*T_0 + 2*T_1 + 3*T_2``
+    while [[1,2],[1,2]] represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) +
+    2*T_0(x)*T_1(y) + 2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is
+    ``y``.
 
     Parameters
     ----------
-    cs: array_like
-        1-d array of Chebyshev series coefficients ordered from low to high.
+    c: array_like
+        Array of Chebyshev series coefficients. If c is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     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)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -894,55 +901,78 @@ def chebder(cs, m=1, scl=1) :
     Examples
     --------
     >>> from numpy.polynomial import chebyshev as C
-    >>> cs = (1,2,3,4)
-    >>> C.chebder(cs)
+    >>> c = (1,2,3,4)
+    >>> C.chebder(c)
     array([ 14.,  12.,  24.])
-    >>> C.chebder(cs,3)
+    >>> C.chebder(c,3)
     array([ 96.])
-    >>> C.chebder(cs,scl=-1)
+    >>> C.chebder(c,scl=-1)
     array([-14., -12., -24.])
-    >>> C.chebder(cs,2,-1)
+    >>> C.chebder(c,2,-1)
     array([ 12.,  96.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of derivation must be integer")
-    if cnt < 0 :
+    if cnt < 0:
         raise ValueError("The order of derivation must be non-negative")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
-    elif cnt >= len(cs):
-        return cs[:1]*0
-    else :
-        zs = _cseries_to_zseries(cs)
-        for i in range(cnt):
-            zs = _zseries_der(zs)*scl
-        return _zseries_to_cseries(zs)
-
+        return c
 
-def chebint(cs, m=1, k=[], lbnd=0, scl=1):
+    c = np.rollaxis(c, iaxis)
+    n = len(c)
+    if cnt >= n:
+        c = c[:1]*0
+    else:
+        for i in range(cnt):
+            n = n - 1
+            c *= scl
+            der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
+            for j in range(n, 2, -1):
+                der[j - 1] = (2*j)*c[j]
+                c[j - 2] += (j*c[j])/(j - 2)
+            if n > 1:
+                der[1] = 4*c[2]
+            der[0] = c[1]
+            c = der
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
+
+
+def chebint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
     """
     Integrate a Chebyshev series.
 
-    Returns, as a C-series, the input C-series `cs`, integrated `m` times
-    from `lbnd` to `x`.  At each iteration the resulting series is
+    Returns the Chebeshev series coefficients `c` integrated `m` times from
+    `lbnd` along `axis`. 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 C-series "term" to highest, e.g.,
-    [1,2,3] represents the series :math:`T_0(x) + 2T_1(x) + 3T_2(x)`.
+    see the Notes section below.)  The argument `c` is an array of
+    coefficients from low to high degree along each axix, e.g., [1,2,3]
+    represents the series ``T_0 + 2*T_1 + 3*T_2`` while [[1,2],[1,2]]
+    represents ``1*T_0(x)*T_0(y) + 1*T_1(x)*T_0(y) + 2*T_0(x)*T_1(y) +
+    2*T_1(x)*T_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
 
     Parameters
     ----------
-    cs : array_like
-        1-d array of C-series coefficients, ordered from low to high.
+    c : array_like
+        Array of Chebyshev series coefficients. If c is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     m : int, optional
         Order of integration, must be positive. (Default: 1)
     k : {[], list, scalar}, optional
@@ -956,6 +986,8 @@ def chebint(cs, m=1, k=[], lbnd=0, scl=1):
     scl : scalar, optional
         Following each integration the result is *multiplied* by `scl`
         before the integration constant is added. (Default: 1)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -988,23 +1020,26 @@ def chebint(cs, m=1, k=[], lbnd=0, scl=1):
     Examples
     --------
     >>> from numpy.polynomial import chebyshev as C
-    >>> cs = (1,2,3)
-    >>> C.chebint(cs)
+    >>> c = (1,2,3)
+    >>> C.chebint(c)
     array([ 0.5, -0.5,  0.5,  0.5])
-    >>> C.chebint(cs,3)
+    >>> C.chebint(c,3)
     array([ 0.03125   , -0.1875    ,  0.04166667, -0.05208333,  0.01041667,
             0.00625   ])
-    >>> C.chebint(cs, k=3)
+    >>> C.chebint(c, k=3)
     array([ 3.5, -0.5,  0.5,  0.5])
-    >>> C.chebint(cs,lbnd=-2)
+    >>> C.chebint(c,lbnd=-2)
     array([ 8.5, -0.5,  0.5,  0.5])
-    >>> C.chebint(cs,scl=-2)
+    >>> C.chebint(c,scl=-2)
     array([-1.,  1., -1., -1.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if not np.iterable(k):
         k = [k]
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of integration must be integer")
@@ -1012,24 +1047,37 @@ def chebint(cs, m=1, k=[], lbnd=0, scl=1):
         raise ValueError("The order of integration must be non-negative")
     if len(k) > cnt :
         raise ValueError("Too many integration constants")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
+        return c
 
+    c = np.rollaxis(c, iaxis)
     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]
+        n = len(c)
+        c *= scl
+        if n == 1 and np.all(c[0] == 0):
+            c[0] += k[i]
         else:
-            zs = _cseries_to_zseries(cs)
-            zs = _zseries_int(zs)
-            cs = _zseries_to_cseries(zs)
-            cs[0] += k[i] - chebval(lbnd, cs)
-    return cs
+            tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
+            tmp[0] = c[0]*0
+            tmp[1] = c[0]
+            if n > 1:
+                tmp[2] = c[1]/4
+            for j in range(2, n):
+                t = c[j]/(2*j + 1)
+                tmp[j + 1] = c[j]/(2*(j + 1))
+                tmp[j - 1] -= c[j]/(2*(j - 1))
+            tmp[0] += k[i] - chebval(lbnd, tmp)
+            c = tmp
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
 def chebval(x, c, tensor=True):
@@ -1094,9 +1142,9 @@ def chebval(x, c, tensor=True):
     --------
 
     """
-    c = np.array(c, ndmin=1, copy=0)
-    if c.dtype.char not in 'efdgFDGO':
-        c = c + 0.0
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if isinstance(x, (tuple, list)):
         x = np.asarray(x)
     if isinstance(x, np.ndarray) and tensor:
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 801e8f8d9..340e6009b 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -625,26 +625,33 @@ def hermpow(cs, pow, maxpower=16) :
         return prd
 
 
-def hermder(cs, m=1, scl=1) :
+def hermder(c, m=1, scl=1, axis=0) :
     """
     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``.
+    Returns the Hermite series coefficients `c` differentiated `m` times
+    along `axis`.  At each iteration the result is multiplied by `scl` (the
+    scaling factor is for use in a linear change of variable). The argument
+    `c` is an array of coefficients from low to high degree along each
+    axis, e.g., [1,2,3] represents the series ``1*H_0 + 2*H_1 + 3*H_2``
+    while [[1,2],[1,2]] represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) +
+    2*H_0(x)*H_1(y) + 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is
+    ``y``.
 
     Parameters
     ----------
-    cs: array_like
-        1-d array of Hermite series coefficients ordered from low to high.
+    c : array_like
+        Array of Hermite series coefficients. If `c` is multidimensional the
+        different axis correspond to different variables with the degree in
+        each axis given by the corresponding index.
     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)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -671,48 +678,63 @@ def hermder(cs, m=1, scl=1) :
     array([ 1.,  2.,  3.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of derivation must be integer")
-    if cnt < 0 :
+    if cnt < 0:
         raise ValueError("The order of derivation must be non-negative")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
-    elif cnt >= len(cs):
-        return cs[:1]*0
+        return c
+
+    c = np.rollaxis(c, iaxis)
+    n = len(c)
+    if cnt >= n:
+        c = c[:1]*0
     else :
         for i in range(cnt):
-            n = len(cs) - 1
-            cs *= scl
-            der = np.empty(n, dtype=cs.dtype)
+            n = n - 1
+            c *= scl
+            der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
             for j in range(n, 0, -1):
-                der[j - 1] = (2*j)*cs[j]
-            cs = der
-        return cs
+                der[j - 1] = (2*j)*c[j]
+            c = der
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
-def hermint(cs, m=1, k=[], lbnd=0, scl=1):
+def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
     """
     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.
+    Returns the Hermite series coefficients `c` integrated `m` times from
+    `lbnd` along `axis`. 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)`.
+    see the Notes section below.)  The argument `c` is an array of
+    coefficients from low to high degree along each axix, e.g., [1,2,3]
+    represents the series ``H_0 + 2*H_1 + 3*H_2`` while [[1,2],[1,2]]
+    represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) + 2*H_0(x)*H_1(y) +
+    2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
 
     Parameters
     ----------
-    cs : array_like
-        1-d array of Hermite series coefficients, ordered from low to high.
+    c : array_like
+        Array of Hermite series coefficients. If c is multidimensional the
+        different axis correspond to different variables with the degree in
+        each axis given by the corresponding index.
     m : int, optional
         Order of integration, must be positive. (Default: 1)
     k : {[], list, scalar}, optional
@@ -726,6 +748,8 @@ def hermint(cs, m=1, k=[], lbnd=0, scl=1):
     scl : scalar, optional
         Following each integration the result is *multiplied* by `scl`
         before the integration constant is added. (Default: 1)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -760,7 +784,7 @@ def hermint(cs, m=1, k=[], lbnd=0, scl=1):
     >>> from numpy.polynomial.hermite import hermint
     >>> hermint([1,2,3]) # integrate once, value 0 at 0.
     array([ 1. ,  0.5,  0.5,  0.5])
-    >>> hermint([1,2,3], m=2) # integrate twice, value & deriv 0 at 0 
+    >>> hermint([1,2,3], m=2) # integrate twice, value & deriv 0 at 0
     array([-0.5       ,  0.5       ,  0.125     ,  0.08333333,  0.0625    ])
     >>> hermint([1,2,3], k=1) # integrate once, value 1 at 0.
     array([ 2. ,  0.5,  0.5,  0.5])
@@ -770,9 +794,12 @@ def hermint(cs, m=1, k=[], lbnd=0, scl=1):
     array([ 1.66666667, -0.5       ,  0.125     ,  0.08333333,  0.0625    ])
 
     """
-    cnt = int(m)
-    if np.isscalar(k) :
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    if not np.iterable(k):
         k = [k]
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of integration must be integer")
@@ -780,27 +807,33 @@ def hermint(cs, m=1, k=[], lbnd=0, scl=1):
         raise ValueError("The order of integration must be non-negative")
     if len(k) > cnt :
         raise ValueError("Too many integration constants")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
+        return c
 
+    c = np.rollaxis(c, iaxis)
     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]
+        n = len(c)
+        c *= scl
+        if n == 1 and np.all(c[0] == 0):
+            c[0] += k[i]
         else:
-            tmp = np.empty(n + 1, dtype=cs.dtype)
-            tmp[0] = cs[0]*0
-            tmp[1] = cs[0]/2
+            tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
+            tmp[0] = c[0]*0
+            tmp[1] = c[0]/2
             for j in range(1, n):
-                tmp[j + 1] = cs[j]/(2*(j + 1))
+                tmp[j + 1] = c[j]/(2*(j + 1))
             tmp[0] += k[i] - hermval(lbnd, tmp)
-            cs = tmp
-    return cs
+            c = tmp
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
 def hermval(x, c, tensor=True):
@@ -873,8 +906,8 @@ def hermval(x, c, tensor=True):
 
     """
     c = np.array(c, ndmin=1, copy=0)
-    if c.dtype.char not in 'efdgFDGO':
-        c = c + 0.0
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if isinstance(x, (tuple, list)):
         x = np.asarray(x)
     if isinstance(x, np.ndarray) and tensor:
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index 8307cf80c..ab9e397be 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -624,26 +624,33 @@ def hermepow(cs, pow, maxpower=16) :
         return prd
 
 
-def hermeder(cs, m=1, scl=1) :
+def hermeder(c, m=1, scl=1, axis=0) :
     """
-    Differentiate a Hermite series.
+    Differentiate a Hermite_e 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``.
+    Returns the series coefficients `c` differentiated `m` times along
+    `axis`.  At each iteration the result is multiplied by `scl` (the
+    scaling factor is for use in a linear change of variable). The argument
+    `c` is an array of coefficients from low to high degree along each
+    axis, e.g., [1,2,3] represents the series ``1*He_0 + 2*He_1 + 3*He_2``
+    while [[1,2],[1,2]] represents ``1*He_0(x)*He_0(y) + 1*He_1(x)*He_0(y)
+    + 2*He_0(x)*He_1(y) + 2*He_1(x)*He_1(y)`` if axis=0 is ``x`` and axis=1
+    is ``y``.
 
     Parameters
     ----------
-    cs: array_like
-        1-d array of Hermite series coefficients ordered from low to high.
+    c: array_like
+        Array of Hermite_e series coefficients. If `c` is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     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)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -670,48 +677,63 @@ def hermeder(cs, m=1, scl=1) :
     array([ 1.,  2.,  3.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of derivation must be integer")
-    if cnt < 0 :
+    if cnt < 0:
         raise ValueError("The order of derivation must be non-negative")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
-    elif cnt >= len(cs):
-        return cs[:1]*0
+        return c
+
+    c = np.rollaxis(c, iaxis)
+    n = len(c)
+    if cnt >= n:
+        return c[:1]*0
     else :
         for i in range(cnt):
-            n = len(cs) - 1
-            cs *= scl
-            der = np.empty(n, dtype=cs.dtype)
+            n = n - 1
+            c *= scl
+            der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
             for j in range(n, 0, -1):
-                der[j - 1] = j*cs[j]
-            cs = der
-        return cs
+                der[j - 1] = j*c[j]
+            c = der
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
-def hermeint(cs, m=1, k=[], lbnd=0, scl=1):
+def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
     """
-    Integrate a Hermite series.
+    Integrate a Hermite_e 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.
+    Returns the Hermite_e series coefficients `c` integrated `m` times from
+    `lbnd` along `axis`. 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)`.
+    see the Notes section below.)  The argument `c` is an array of
+    coefficients from low to high degree along each axix, e.g., [1,2,3]
+    represents the series ``H_0 + 2*H_1 + 3*H_2`` while [[1,2],[1,2]]
+    represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) + 2*H_0(x)*H_1(y) +
+    2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
 
     Parameters
     ----------
-    cs : array_like
-        1-d array of Hermite series coefficients, ordered from low to high.
+    c : array_like
+        Array of Hermite_e series coefficients. If c is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     m : int, optional
         Order of integration, must be positive. (Default: 1)
     k : {[], list, scalar}, optional
@@ -725,11 +747,13 @@ def hermeint(cs, m=1, k=[], lbnd=0, scl=1):
     scl : scalar, optional
         Following each integration the result is *multiplied* by `scl`
         before the integration constant is added. (Default: 1)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
     S : ndarray
-        Hermite series coefficients of the integral.
+        Hermite_e series coefficients of the integral.
 
     Raises
     ------
@@ -765,13 +789,16 @@ def hermeint(cs, m=1, k=[], lbnd=0, scl=1):
     array([ 2.,  1.,  1.,  1.])
     >>> hermeint([1, 2, 3], lbnd=-1) # integrate once, value 0 at -1
     array([-1.,  1.,  1.,  1.])
-    >>> hermeint([1, 2, 3], m=2, k=[1,2], lbnd=-1)
+    >>> hermeint([1, 2, 3], m=2, k=[1, 2], lbnd=-1)
     array([ 1.83333333,  0.        ,  0.5       ,  0.33333333,  0.25      ])
 
     """
-    cnt = int(m)
-    if np.isscalar(k) :
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    if not np.iterable(k):
         k = [k]
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of integration must be integer")
@@ -779,27 +806,33 @@ def hermeint(cs, m=1, k=[], lbnd=0, scl=1):
         raise ValueError("The order of integration must be non-negative")
     if len(k) > cnt :
         raise ValueError("Too many integration constants")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
+        return c
 
+    c = np.rollaxis(c, iaxis)
     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]
+        n = len(c)
+        c *= scl
+        if n == 1 and np.all(c[0] == 0):
+            c[0] += k[i]
         else:
-            tmp = np.empty(n + 1, dtype=cs.dtype)
-            tmp[0] = cs[0]*0
-            tmp[1] = cs[0]
+            tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
+            tmp[0] = c[0]*0
+            tmp[1] = c[0]
             for j in range(1, n):
-                tmp[j + 1] = cs[j]/(j + 1)
+                tmp[j + 1] = c[j]/(j + 1)
             tmp[0] += k[i] - hermeval(lbnd, tmp)
-            cs = tmp
-    return cs
+            c = tmp
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
 def hermeval(x, c, tensor=True):
@@ -871,8 +904,8 @@ def hermeval(x, c, tensor=True):
 
     """
     c = np.array(c, ndmin=1, copy=0)
-    if c.dtype.char not in 'efdgFDGO':
-        c = c + 0.0
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if isinstance(x, (tuple, list)):
         x = np.asarray(x)
     if isinstance(x, np.ndarray) and tensor:
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 39463b0cb..07f943bb1 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -621,26 +621,33 @@ def lagpow(cs, pow, maxpower=16) :
         return prd
 
 
-def lagder(cs, m=1, scl=1) :
+def lagder(c, m=1, scl=1, axis=0) :
     """
     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``.
+    Returns the Laguerre series coefficients `c` differentiated `m` times
+    along `axis`.  At each iteration the result is multiplied by `scl` (the
+    scaling factor is for use in a linear change of variable). The argument
+    `c` is an array of coefficients from low to high degree along each
+    axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2``
+    while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) +
+    2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is
+    ``y``.
 
     Parameters
     ----------
-    cs: array_like
-        1-d array of Laguerre series coefficients ordered from low to high.
+    c: array_like
+        Array of Laguerre series coefficients. If `c` is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     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)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -667,49 +674,66 @@ def lagder(cs, m=1, scl=1) :
     array([ 1.,  2.,  3.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of derivation must be integer")
-    if cnt < 0 :
+    if cnt < 0:
         raise ValueError("The order of derivation must be non-negative")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
-    elif cnt >= len(cs):
-        return cs[:1]*0
+        return c
+
+    c = np.rollaxis(c, iaxis)
+    n = len(c)
+    if cnt >= n:
+        c = c[: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):
+            n = n - 1
+            c *= scl
+            der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
+            for j in range(n, 1, -1):
+                der[j - 1] = -c[j]
+                c[j - 1] += c[j]
+            der[0] = -c[1]
+            c = der
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
+
+
+def lagint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
     """
     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.
+    Returns the Laguerre series coefficients `c` integrated `m` times from
+    `lbnd` along `axis`. 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)`.
+    see the Notes section below.)  The argument `c` is an array of
+    coefficients from low to high degree along each axix, e.g., [1,2,3]
+    represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]]
+    represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) +
+    2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
+
 
     Parameters
     ----------
-    cs : array_like
-        1-d array of Laguerre series coefficients, ordered from low to high.
+    c : array_like
+        Array of Laguerre series coefficients. If `c` is multidimensional
+        the different axis correspond to different variables with the
+        degree in each axis given by the corresponding index.
     m : int, optional
         Order of integration, must be positive. (Default: 1)
     k : {[], list, scalar}, optional
@@ -723,6 +747,8 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1):
     scl : scalar, optional
         Following each integration the result is *multiplied* by `scl`
         before the integration constant is added. (Default: 1)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -767,9 +793,12 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1):
     array([ 11.16666667,  -5.        ,  -3.        ,   2.        ])
 
     """
-    cnt = int(m)
-    if np.isscalar(k) :
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    if not np.iterable(k):
         k = [k]
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of integration must be integer")
@@ -777,28 +806,34 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1):
         raise ValueError("The order of integration must be non-negative")
     if len(k) > cnt :
         raise ValueError("Too many integration constants")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
+        return c
 
+    c = np.rollaxis(c, iaxis)
     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]
+        n = len(c)
+        c *= scl
+        if n == 1 and np.all(c[0] == 0):
+            c[0] += k[i]
         else:
-            tmp = np.empty(n + 1, dtype=cs.dtype)
-            tmp[0] = cs[0]
-            tmp[1] = -cs[0]
+            tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
+            tmp[0] = c[0]
+            tmp[1] = -c[0]
             for j in range(1, n):
-                tmp[j] += cs[j]
-                tmp[j + 1] = -cs[j]
+                tmp[j] += c[j]
+                tmp[j + 1] = -c[j]
             tmp[0] += k[i] - lagval(lbnd, tmp)
-            cs = tmp
-    return cs
+            c = tmp
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
 def lagval(x, c, tensor=True):
@@ -871,8 +906,8 @@ def lagval(x, c, tensor=True):
 
     """
     c = np.array(c, ndmin=1, copy=0)
-    if c.dtype.char not in 'efdgFDGO':
-        c = c + 0.0
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if isinstance(x, (tuple, list)):
         x = np.asarray(x)
     if isinstance(x, np.ndarray) and tensor:
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 00dbacebe..8ef439ab4 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -632,26 +632,33 @@ def legpow(cs, pow, maxpower=16) :
         return prd
 
 
-def legder(cs, m=1, scl=1) :
+def legder(c, m=1, scl=1, axis=0) :
     """
     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``.
+    Returns the Legendre series coefficients `c` differentiated `m` times
+    along `axis`.  At each iteration the result is multiplied by `scl` (the
+    scaling factor is for use in a linear change of variable). The argument
+    `c` is an array of coefficients from low to high degree along each
+    axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2``
+    while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) +
+    2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is
+    ``y``.
 
     Parameters
     ----------
-    cs : array_like
-        1-D array of Legendre series coefficients ordered from low to high.
+    c : array_like
+        Array of Legendre series coefficients. If c is multidimensional the
+        different axis correspond to different variables with the degree in
+        each axis given by the corresponding index.
     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)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
@@ -672,60 +679,78 @@ def legder(cs, m=1, scl=1) :
     Examples
     --------
     >>> from numpy.polynomial import legendre as L
-    >>> cs = (1,2,3,4)
-    >>> L.legder(cs)
+    >>> c = (1,2,3,4)
+    >>> L.legder(c)
     array([  6.,   9.,  20.])
-    >>> L.legder(cs,3)
+    >>> L.legder(c, 3)
     array([ 60.])
-    >>> L.legder(cs,scl=-1)
+    >>> L.legder(c, scl=-1)
     array([ -6.,  -9., -20.])
-    >>> L.legder(cs,2,-1)
+    >>> L.legder(c, 2,-1)
     array([  9.,  60.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of derivation must be integer")
-    if cnt < 0 :
+    if cnt < 0:
         raise ValueError("The order of derivation must be non-negative")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
-    elif cnt >= len(cs):
-        return cs[:1]*0
+        return c
+
+    c = np.rollaxis(c, iaxis)
+    n = len(c)
+    if cnt >= n:
+        c = c[: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):
+            n = n - 1
+            c *= scl
+            der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
+            for j in range(n, 2, -1):
+                der[j - 1] = (2*j - 1)*c[j]
+                c[j - 2] += c[j]
+            if n > 1:
+                der[1] = 3*c[2]
+            der[0] = c[1]
+            c = der
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
+
+
+def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
     """
     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.
+    Returns the Legendre series coefficients `c` integrated `m` times from
+    `lbnd` along `axis`. 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)`.
+    see the Notes section below.)  The argument `c` is an array of
+    coefficients from low to high degree along each axix, e.g., [1,2,3]
+    represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]]
+    represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) +
+    2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
 
     Parameters
     ----------
-    cs : array_like
-        1-d array of Legendre series coefficients, ordered from low to high.
+    c : array_like
+        Array of Legendre series coefficients. If c is multidimensional the
+        different axis correspond to different variables with the degree in
+        each axis given by the corresponding index.
     m : int, optional
         Order of integration, must be positive. (Default: 1)
     k : {[], list, scalar}, optional
@@ -739,11 +764,13 @@ def legint(cs, m=1, k=[], lbnd=0, scl=1):
     scl : scalar, optional
         Following each integration the result is *multiplied* by `scl`
         before the integration constant is added. (Default: 1)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
     S : ndarray
-        Legendre series coefficients of the integral.
+        Legendre series coefficient array of the integral.
 
     Raises
     ------
@@ -771,23 +798,26 @@ def legint(cs, m=1, k=[], lbnd=0, scl=1):
     Examples
     --------
     >>> from numpy.polynomial import legendre as L
-    >>> cs = (1,2,3)
+    >>> c = (1,2,3)
     >>> L.legint(cs)
     array([ 0.33333333,  0.4       ,  0.66666667,  0.6       ])
-    >>> L.legint(cs,3)
+    >>> L.legint(c, 3)
     array([  1.66666667e-02,  -1.78571429e-02,   4.76190476e-02,
             -1.73472348e-18,   1.90476190e-02,   9.52380952e-03])
-    >>> L.legint(cs, k=3)
+    >>> L.legint(c, k=3)
     array([ 3.33333333,  0.4       ,  0.66666667,  0.6       ])
-    >>> L.legint(cs, lbnd=-2)
+    >>> L.legint(c, lbnd=-2)
     array([ 7.33333333,  0.4       ,  0.66666667,  0.6       ])
-    >>> L.legint(cs, scl=2)
+    >>> L.legint(c, scl=2)
     array([ 0.66666667,  0.8       ,  1.33333333,  1.2       ])
 
     """
-    cnt = int(m)
-    if np.isscalar(k) :
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    if not np.iterable(k):
         k = [k]
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of integration must be integer")
@@ -795,29 +825,37 @@ def legint(cs, m=1, k=[], lbnd=0, scl=1):
         raise ValueError("The order of integration must be non-negative")
     if len(k) > cnt :
         raise ValueError("Too many integration constants")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
+        return c
 
+    c = np.rollaxis(c, iaxis)
     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]
+        n = len(c)
+        c *= scl
+        if n == 1 and np.all(c[0] == 0):
+            c[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 = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
+            tmp[0] = c[0]*0
+            tmp[1] = c[0]
+            if n > 1:
+                tmp[2] = c[1]/3
+            for j in range(2, n):
+                t = c[j]/(2*j + 1)
                 tmp[j + 1] = t
                 tmp[j - 1] -= t
             tmp[0] += k[i] - legval(lbnd, tmp)
-            cs = tmp
-    return cs
+            c = tmp
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
 def legval(x, c, tensor=True):
@@ -882,12 +920,12 @@ def legval(x, c, tensor=True):
 
     """
     c = np.array(c, ndmin=1, copy=0)
-    if c.dtype.char not in 'efdgFDGO':
-        c = c + 0.0
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if isinstance(x, (tuple, list)):
         x = np.asarray(x)
     if isinstance(x, np.ndarray) and tensor:
-       c = c.reshape(c.shape + (1,)*x.ndim)
+        c = c.reshape(c.shape + (1,)*x.ndim)
 
     if len(c) == 1 :
         c0 = c[0]
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index bf9c1a216..87acfa616 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -414,8 +414,8 @@ def polypow(cs, pow, maxpower=None) :
     Parameters
     ----------
     cs : array_like
-        1d array of chebyshev series coefficients ordered from low to
-        high.
+        1d array of array of series coefficients ordered from low to
+        high degree.
     pow : integer
         Power to which the series will be raised
     maxpower : integer, optional
@@ -425,11 +425,11 @@ def polypow(cs, pow, maxpower=None) :
     Returns
     -------
     coef : ndarray
-        Chebyshev series of power.
+        Power series of power.
 
     See Also
     --------
-    chebadd, chebsub, chebmul, chebdiv
+    polyadd, polysub, polymul, polydiv
 
     Examples
     --------
@@ -455,31 +455,37 @@ def polypow(cs, pow, maxpower=None) :
         return prd
 
 
-def polyder(cs, m=1, scl=1):
+def polyder(c, m=1, scl=1, axis=0):
     """
     Differentiate a polynomial.
 
-    Returns the polynomial `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 polynomial ``1 + 2*x + 3*x**2``.
+    Returns the polynomial coefficients `c` differentiated `m` times along
+    `axis`.  At each iteration the result is multiplied by `scl` (the
+    scaling factor is for use in a linear change of variable).  The
+    argument `c` is an array of coefficients from low to high degree along
+    each axis, e.g., [1,2,3] represents the polynomial ``1 + 2*x + 3*x**2``
+    while [[1,2],[1,2]] represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is
+    ``x`` and axis=1 is ``y``.
 
     Parameters
     ----------
-    cs: array_like
-        1-d array of polynomial coefficients ordered from low to high.
+    c: array_like
+        Array of polynomial coefficients. If c is multidimensional the
+        different axis correspond to different variables with the degree
+        in each axis given by the corresponding index.
     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)
+    axis : int, optional
+        Axis over which the derivative is taken. (Default: 0).
 
     Returns
     -------
     der : ndarray
-        Polynomial of the derivative.
+        Polynomial coefficients of the derivative.
 
     See Also
     --------
@@ -488,55 +494,71 @@ def polyder(cs, m=1, scl=1):
     Examples
     --------
     >>> from numpy import polynomial as P
-    >>> cs = (1,2,3,4) # 1 + 2x + 3x**2 + 4x**3
-    >>> P.polyder(cs) # (d/dx)(cs) = 2 + 6x + 12x**2
+    >>> c = (1,2,3,4) # 1 + 2x + 3x**2 + 4x**3
+    >>> P.polyder(c) # (d/dx)(c) = 2 + 6x + 12x**2
     array([  2.,   6.,  12.])
-    >>> P.polyder(cs,3) # (d**3/dx**3)(cs) = 24
+    >>> P.polyder(c,3) # (d**3/dx**3)(c) = 24
     array([ 24.])
-    >>> P.polyder(cs,scl=-1) # (d/d(-x))(cs) = -2 - 6x - 12x**2
+    >>> P.polyder(c,scl=-1) # (d/d(-x))(c) = -2 - 6x - 12x**2
     array([ -2.,  -6., -12.])
-    >>> P.polyder(cs,2,-1) # (d**2/d(-x)**2)(cs) = 6 + 24x
+    >>> P.polyder(c,2,-1) # (d**2/d(-x)**2)(c) = 6 + 24x
     array([  6.,  24.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     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")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
-    elif cnt >= len(cs):
-        return cs[:1]*0
+        return c
+
+    c = np.rollaxis(c, iaxis)
+    n = len(c)
+    if cnt >= n:
+        c = c[:1]*0
     else :
-        n = len(cs)
-        d = np.arange(n)*scl
         for i in range(cnt):
-            cs[i:] *= d[:n-i]
-        return cs[i+1:].copy()
+            n = n - 1
+            c *= scl
+            der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
+            for j in range(n, 0, -1):
+                der[j - 1] = j*c[j]
+            c = der
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
-def polyint(cs, m=1, k=[], lbnd=0, scl=1):
+def polyint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
     """
     Integrate a polynomial.
 
-    Returns the polynomial `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
-    term to highest, e.g., [1,2,3] represents the polynomial
-    ``1 + 2*x + 3*x**2``.
+    Returns the polynomial coefficients `c` integrated `m` times from
+    `lbnd` along `axis`.  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 `c` is an array of
+    coefficients, from low to high degree along each axis, e.g., [1,2,3]
+    represents the polynomial ``1 + 2*x + 3*x**2`` while [[1,2],[1,2]]
+    represents ``1 + 1*x + 2*y + 2*x*y`` if axis=0 is ``x`` and axis=1 is
+    ``y``.
 
     Parameters
     ----------
-    cs : array_like
+    c : array_like
         1-d array of polynomial coefficients, ordered from low to high.
     m : int, optional
         Order of integration, must be positive. (Default: 1)
@@ -551,11 +573,13 @@ def polyint(cs, m=1, k=[], lbnd=0, scl=1):
     scl : scalar, optional
         Following each integration the result is *multiplied* by `scl`
         before the integration constant is added. (Default: 1)
+    axis : int, optional
+        Axis over which the integral is taken. (Default: 0).
 
     Returns
     -------
     S : ndarray
-        Coefficients of the integral.
+        Coefficient array of the integral.
 
     Raises
     ------
@@ -577,23 +601,26 @@ def polyint(cs, m=1, k=[], lbnd=0, scl=1):
     Examples
     --------
     >>> from numpy import polynomial as P
-    >>> cs = (1,2,3)
-    >>> P.polyint(cs) # should return array([0, 1, 1, 1])
+    >>> c = (1,2,3)
+    >>> P.polyint(c) # should return array([0, 1, 1, 1])
     array([ 0.,  1.,  1.,  1.])
-    >>> P.polyint(cs,3) # should return array([0, 0, 0, 1/6, 1/12, 1/20])
+    >>> P.polyint(c,3) # should return array([0, 0, 0, 1/6, 1/12, 1/20])
     array([ 0.        ,  0.        ,  0.        ,  0.16666667,  0.08333333,
             0.05      ])
-    >>> P.polyint(cs,k=3) # should return array([3, 1, 1, 1])
+    >>> P.polyint(c,k=3) # should return array([3, 1, 1, 1])
     array([ 3.,  1.,  1.,  1.])
-    >>> P.polyint(cs,lbnd=-2) # should return array([6, 1, 1, 1])
+    >>> P.polyint(c,lbnd=-2) # should return array([6, 1, 1, 1])
     array([ 6.,  1.,  1.,  1.])
-    >>> P.polyint(cs,scl=-2) # should return array([0, -2, -2, -2])
+    >>> P.polyint(c,scl=-2) # should return array([0, -2, -2, -2])
     array([ 0., -2., -2., -2.])
 
     """
-    cnt = int(m)
+    c = np.array(c, ndmin=1, copy=1)
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if not np.iterable(k):
         k = [k]
+    cnt, iaxis = [int(t) for t in [m, axis]]
 
     if cnt != m:
         raise ValueError("The order of integration must be integer")
@@ -601,25 +628,34 @@ def polyint(cs, m=1, k=[], lbnd=0, scl=1):
         raise ValueError("The order of integration must be non-negative")
     if len(k) > cnt :
         raise ValueError("Too many integration constants")
+    if iaxis != axis:
+        raise ValueError("The axis must be integer")
+    if not -c.ndim <= iaxis < c.ndim:
+        raise ValueError("The axis is out of range")
+    if iaxis < 0:
+        iaxis += c.ndim
+
 
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
     if cnt == 0:
-        return cs
+        return c
 
     k = list(k) + [0]*(cnt - len(k))
+    c = np.rollaxis(c, iaxis)
     for i in range(cnt):
-        n = len(cs)
-        cs *= scl
-        if n == 1 and cs[0] == 0:
-            cs[0] += k[i]
+        n = len(c)
+        c *= scl
+        if n == 1 and np.all(c[0] == 0):
+            c[0] += k[i]
         else:
-            tmp = np.empty(n + 1, dtype=cs.dtype)
-            tmp[0] = cs[0]*0
-            tmp[1:] = cs/np.arange(1, n + 1)
+            tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
+            tmp[0] = c[0]*0
+            tmp[1] = c[0]
+            for j in range(1, n):
+                tmp[j + 1] = c[j]/(j + 1)
             tmp[0] += k[i] - polyval(lbnd, tmp)
-            cs = tmp
-    return cs
+            c = tmp
+    c = np.rollaxis(c, 0, iaxis + 1)
+    return c
 
 
 def polyval(x, c, tensor=True):
@@ -704,12 +740,12 @@ def polyval(x, c, tensor=True):
 
     """
     c = np.array(c, ndmin=1, copy=0)
-    if c.dtype.char not in 'efdgFDGO':
-        c = c + 0.0
+    if c.dtype.char in '?bBhHiIlLqQpP':
+        c = c.astype(np.double)
     if isinstance(x, (tuple, list)):
         x = np.asarray(x)
     if isinstance(x, np.ndarray) and tensor:
-       c = c.reshape(c.shape + (1,)*x.ndim)
+        c = c.reshape(c.shape + (1,)*x.ndim)
 
     c0 = c[-1] + x*0
     for i in range(2, len(c) + 1) :
-- 
cgit v1.2.1