ENH: Modify the various polynomial series so that multidimensional coefficient arrays can be used. Add functions for evaluation of 2D and 3D polynomial series evaluated either on a specified set of points or on a cartesian product of 1D points.

The new functions have names polyval2d, polygrid2d, polyval3d, and
polygrid3d, where 'poly' can be replaced by any of 'leg', 'cheb', 'lag',
'herm', or 'herme'. These additional functions should cover the common
multidimensional cases and provide examples for anyone who wants to go to
higher dimensions.

6 files changed, 1361 insertions, 222 deletions

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index b4f50d90e..6212f2bc5 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -25,6 +25,10 @@ Arithmetic
 - `chebdiv` -- divide one Chebyshev series by another.
 - `chebpow` -- raise a Chebyshev series to an positive integer power
 - `chebval` -- evaluate a Chebyshev series at given points.
+- `chebval2d` -- evaluate a 2D Chebyshev series at given points.
+- `chebval3d` -- evaluate a 3D Chebyshev series at given points.
+- `chebgrid2d` -- evaluate a 2D Chebyshev series on a Cartesian product.
+- `chebgrid3d` -- evaluate a 3D Chebyshev series on a Cartesian product.
 
 Calculus
 --------
@@ -78,18 +82,19 @@ References
 """
 from __future__ import division
 
-__all__ = ['chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline',
-        'chebadd', 'chebsub', 'chebmulx', 'chebmul', 'chebdiv', 'chebpow',
-        'chebval', 'chebder', 'chebint', 'cheb2poly', 'poly2cheb',
-        'chebfromroots', 'chebvander', 'chebfit', 'chebtrim', 'chebroots',
-        'chebpts1', 'chebpts2', 'Chebyshev']
-
 import numpy as np
 import numpy.linalg as la
 import polyutils as pu
 import warnings
 from polytemplate import polytemplate
 
+__all__ = ['chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline',
+    'chebadd', 'chebsub', 'chebmulx', 'chebmul', 'chebdiv', 'chebpow',
+    'chebval', 'chebval2d', 'chebval3d', 'chebgrid2d', 'chebgrid3d',
+    'chebder', 'chebint', 'cheb2poly', 'poly2cheb', 'chebfromroots',
+    'chebvander', 'chebfit', 'chebtrim', 'chebroots', 'chebpts1',
+    'chebpts2', 'Chebyshev']
+
 chebtrim = pu.trimcoef
 
 #
@@ -1025,36 +1030,59 @@ def chebint(cs, m=1, k=[], lbnd=0, scl=1):
     return cs
 
 
-def chebval(x, cs):
-    """Evaluate a Chebyshev series.
+def chebval(x, c, tensor=True):
+    """
+    Evaluate a Chebyshev series.
+
+    If `c` is of length ``n + 1``, this function returns the value:
 
-    If `cs` is of length `n`, this function returns :
+    ``p(x) = c[0]*T_0(x) + c[1]*T_1(x) + ... + c[n]*T_n(x)``
 
-    ``p(x) = cs[0]*T_0(x) + cs[1]*T_1(x) + ... + cs[n-1]*T_{n-1}(x)``
+    If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)``
+    will have the same shape as `x`.  If `x` is a algebra_like object that
+    supports multiplication and addition with itself and the values in `c`,
+    then an object of the same type is returned.
 
-    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.
+    In the case where c is multidimensional, the shape of the result
+    depends on the value of `tensor`. If tensor is true the shape of the
+    return will be ``c.shape[1:] + x.shape``, where the shape of a scalar
+    is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if
+    `x` is broadcast compatible with that.
+
+    If there are trailing zeros in the coefficients they still take part in
+    the evaluation, so they should be avoided if efficiency is a concern.
 
     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.
+    x : array_like, algebra_like
+        If x is a list or tuple, it is converted to an ndarray. Otherwise
+        it is left unchanged and if it isn't an ndarray it is treated as a
+        scalar. In either case, `x` or any element of an ndarray must
+        support addition and multiplication with itself and the elements of
+        `c`.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree n are contained in ``c[n]``. If `c` is multidimesional the
+        remaining indices enumerate multiple Chebyshev series. In the two
+        dimensional case the coefficients may be thought of as stored in
+        the columns of `c`.
+    tensor : boolean, optional
+        If true, the coefficient array shape is extended with ones on the
+        right, one for each dimension of `x`. Scalars are treated as having
+        dimension 0 for this action. The effect is that every column of
+        coefficients in `c` is evaluated for every value in `x`. If False,
+        the `x` are broadcast over the columns of `c` in the usual way.
+        This gives some flexibility in evaluations in the multidimensional
+        case. The default value it ``True``.
 
     Returns
     -------
-    values : ndarray, ring_like
-        If the return is an ndarray then it has the same shape as `x`.
+    values : ndarray, algebra_like
+        The shape of the return value is described above.
 
     See Also
     --------
-    chebfit
-
-    Examples
-    --------
+    chebval2d, chebgrid2d, chebval3d, chebgrid3d
 
     Notes
     -----
@@ -1064,28 +1092,186 @@ def chebval(x, cs):
     --------
 
     """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    if isinstance(x, tuple) or isinstance(x, list) :
+    c = np.array(c, ndmin=1, copy=0)
+    if c.dtype.char not in 'efdgFDGO':
+        c = c + 0.0
+    if isinstance(x, (tuple, list)):
         x = np.asarray(x)
+    if isinstance(x, np.ndarray) and tensor:
+       c = c.reshape(c.shape + (1,)*x.ndim)
 
-    if len(cs) == 1 :
-        c0 = cs[0]
+    if len(c) == 1 :
+        c0 = c[0]
         c1 = 0
-    elif len(cs) == 2 :
-        c0 = cs[0]
-        c1 = cs[1]
+    elif len(c) == 2 :
+        c0 = c[0]
+        c1 = c[1]
     else :
         x2 = 2*x
-        c0 = cs[-2]
-        c1 = cs[-1]
-        for i in range(3, len(cs) + 1) :
+        c0 = c[-2]
+        c1 = c[-1]
+        for i in range(3, len(c) + 1) :
             tmp = c0
-            c0 = cs[-i] - c1
+            c0 = c[-i] - c1
             c1 = tmp + c1*x2
     return c0 + c1*x
 
 
+def chebval2d(x, y, c):
+    """
+    Evaluate 2D Chebyshev series at points (x,y).
+
+    This function returns the values:
+
+    ``p(x,y) = \sum_{i,j} c[i,j] * T_i(x) * T_j(y)``
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Chebyshev seres is evaluated at the points
+        ``(x,y)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray.
+        Otherwise it is left unchanged and if it isn't an ndarray it is
+        treated as a scalar.  See `chebval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Chebyshev series at points formed
+        from pairs of corresponding values from `x` and `y`.
+
+    See Also
+    --------
+    chebval, chebgrid2d, chebval3d, chebgrid3d
+
+    """
+    return chebval(y, chebval(x, c), False)
+
+
+def chebgrid2d(x, y, c):
+    """
+    Evaluate 2D Chebyshev series on the Cartesion product of x,y.
+
+    This function returns the values:
+
+    ``p(a,b) = \sum_{i,j} c[i,j] * T_i(a) * T_j(b)``
+
+    where the points ``(a,b)`` consist of all pairs of points formed by
+    taking ``a`` from `x` and ``b`` from `y`. The resulting points form a
+    grid with `x` in the first dimension and `y` in the second.
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Chebyshev series is evaluated at the points in
+        the Cartesian product of `x` and `y`.  If `x` or `y` is a list or
+        tuple, it is first converted to an ndarray, Otherwise it is left
+        unchanged and if it isn't an ndarray it is treated as a scalar. See
+        `chebval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Chebyshev series at points in the
+        Cartesion product of `x` and `y`.
+
+    See Also
+    --------
+    chebval, chebval2d, chebval3d, chebgrid3d
+
+    """
+    return chebval(y, chebval(x, c))
+
+
+def chebval3d(x, y, z, c):
+    """
+    Evaluate 3D Chebyshev series at points (x,y,z).
+
+    This function returns the values:
+
+    ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * T_i(x) * T_j(y) * T_k(z)``
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Chebyshev seres is evaluated at the points
+        ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape.  If
+        any of `x`, `y`, or `z` is a list or tuple, it is first converted
+        to an ndarray. Otherwise it is left unchanged and if it isn't an
+        ndarray it is treated as a scalar.  See `chebval` for explanation of
+        algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Chebyshev series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    chebval, chebval2d, chebgrid2d, chebgrid3d
+
+    """
+    return chebval(z, chebval2d(x, y, c), False)
+
+
+def chebgrid3d(x, y, z, c):
+    """
+    Evaluate 3D Chebyshev series on the Cartesian product of x,y,z.
+
+    This function returns the values:
+
+    ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * T_i(a) * T_j(b) * T_k(c)``
+
+    where the points ``(a,b,c)`` consist of all triples formed by taking
+    ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting
+    points form a grid with `x` in the first dimension, `y` in the second,
+    and `z` in the third.
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Chebyshev seres is evaluated at the points
+        in the cartesian product of `x`, `y`, and `z`
+        ``(x,y,z)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and treated as a scalar.  See
+        `chebval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Chebyshev series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    chebval, chebval2d, chebgrid2d, chebval3d
+
+    """
+    return chebval(z, chebgrid2d(x, y, c))
+
+
 def chebvander(x, deg) :
     """Vandermonde matrix of given degree.
 
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 68713d83a..2fd28a3ff 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -22,6 +22,10 @@ Arithmetic
 - `hermmul` -- multiply two Hermite series.
 - `hermdiv` -- divide one Hermite series by another.
 - `hermval` -- evaluate a Hermite series at given points.
+- `hermval2d` -- evaluate a 2D Hermite series at given points.
+- `hermval3d` -- evaluate a 3D Hermite series at given points.
+- `hermgrid2d` -- evaluate a 2D Hermite series on a Cartesian product.
+- `hermgrid3d` -- evaluate a 3D Hermite series on a Cartesian product.
 
 Calculus
 --------
@@ -50,17 +54,18 @@ See also
 """
 from __future__ import division
 
-__all__ = ['hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline',
-        'hermadd', 'hermsub', 'hermmulx', 'hermmul', 'hermdiv', 'hermval',
-        'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots',
-        'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite']
-
 import numpy as np
 import numpy.linalg as la
 import polyutils as pu
 import warnings
 from polytemplate import polytemplate
 
+__all__ = ['hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline',
+    'hermadd', 'hermsub', 'hermmulx', 'hermmul', 'hermdiv', 'hermpow',
+    'hermval', 'hermval2d', 'hermval3d', 'hermgrid2d', 'hermgrid3d',
+    'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots',
+    'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite']
+
 hermtrim = pu.trimcoef
 
 
@@ -795,36 +800,59 @@ def hermint(cs, m=1, k=[], lbnd=0, scl=1):
     return cs
 
 
-def hermval(x, cs):
-    """Evaluate a Hermite series.
+def hermval(x, c, tensor=True):
+    """
+    Evaluate a Hermite series.
+
+    If `c` is of length ``n + 1``, this function returns the value:
 
-    If `cs` is of length `n`, this function returns :
+    ``p(x) = c[0]*H_0(x) + c[1]*H_1(x) + ... + c[n]*H_n(x)``
 
-    ``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 and `c` is 1 dimensional, then ``p(x)``
+    will have the same shape as `x`.  If `x` is a algebra_like object that
+    supports multiplication and addition with itself and the values in `c`,
+    then an object of the same type is returned.
 
-    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.
+    In the case where c is multidimensional, the shape of the result
+    depends on the value of `tensor`. If tensor is true the shape of the
+    return will be ``c.shape[1:] + x.shape``, where the shape of a scalar
+    is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if
+    `x` is broadcast compatible with that.
+
+    If there are trailing zeros in the coefficients they still take part in
+    the evaluation, so they should be avoided if efficiency is a concern.
 
     Parameters
     ----------
-    x : array_like, ring_like
-        Array of numbers or objects that support multiplication and
-        addition with themselves and with the elements of `cs`.
-    cs : array_like
-        1-d array of Hermite coefficients ordered from low to high.
+    x : array_like, algebra_like
+        If x is a list or tuple, it is converted to an ndarray. Otherwise
+        it is left unchanged and if it isn't an ndarray it is treated as a
+        scalar. In either case, `x` or any element of an ndarray must
+        support addition and multiplication with itself and the elements of
+        `c`.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree n are contained in ``c[n]``. If `c` is multidimesional the
+        remaining indices enumerate multiple Hermite series. In the two
+        dimensional case the coefficients may be thought of as stored in
+        the columns of `c`.
+    tensor : boolean, optional
+        If true, the coefficient array shape is extended with ones on the
+        right, one for each dimension of `x`. Scalars are treated as having
+        dimension 0 for this action. The effect is that every column of
+        coefficients in `c` is evaluated for every value in `x`. If False,
+        the `x` are broadcast over the columns of `c` in the usual way.
+        This gives some flexibility in evaluations in the multidimensional
+        case. The default value it ``True``.
 
     Returns
     -------
-    values : ndarray, ring_like
-        If the return is an ndarray then it has the same shape as `x`.
+    values : ndarray, algebra_like
+        The shape of the return value is described above.
 
     See Also
     --------
-    hermfit
-
-    Examples
-    --------
+    hermval2d, hermgrid2d, hermval3d, hermgrid3d
 
     Notes
     -----
@@ -841,30 +869,188 @@ def hermval(x, cs):
            [ 115.,  203.]])
 
     """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    if isinstance(x, tuple) or isinstance(x, list) :
+    c = np.array(c, ndmin=1, copy=0)
+    if c.dtype.char not in 'efdgFDGO':
+        c = c + 0.0
+    if isinstance(x, (tuple, list)):
         x = np.asarray(x)
+    if isinstance(x, np.ndarray) and tensor:
+       c = c.reshape(c.shape + (1,)*x.ndim)
 
     x2 = x*2
-    if len(cs) == 1 :
-        c0 = cs[0]
+    if len(c) == 1 :
+        c0 = c[0]
         c1 = 0
-    elif len(cs) == 2 :
-        c0 = cs[0]
-        c1 = cs[1]
+    elif len(c) == 2 :
+        c0 = c[0]
+        c1 = c[1]
     else :
-        nd = len(cs)
-        c0 = cs[-2]
-        c1 = cs[-1]
-        for i in range(3, len(cs) + 1) :
+        nd = len(c)
+        c0 = c[-2]
+        c1 = c[-1]
+        for i in range(3, len(c) + 1) :
             tmp = c0
             nd =  nd - 1
-            c0 = cs[-i] - c1*(2*(nd - 1))
+            c0 = c[-i] - c1*(2*(nd - 1))
             c1 = tmp + c1*x2
     return c0 + c1*x2
 
 
+def hermval2d(x, y, c):
+    """
+    Evaluate 2D Hermite series at points (x,y).
+
+    This function returns the values:
+
+    ``p(x,y) = \sum_{i,j} c[i,j] * H_i(x) * H_j(y)``
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Hermite seres is evaluated at the points
+        ``(x,y)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray.
+        Otherwise it is left unchanged and if it isn't an ndarray it is
+        treated as a scalar.  See `hermval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Hermite series at points formed
+        from pairs of corresponding values from `x` and `y`.
+
+    See Also
+    --------
+    hermval, hermgrid2d, hermval3d, hermgrid3d
+
+    """
+    return hermval(y, hermval(x, c), False)
+
+
+def hermgrid2d(x, y, c):
+    """
+    Evaluate 2D Hermite series on the Cartesion product of x,y.
+
+    This function returns the values:
+
+    ``p(a,b) = \sum_{i,j} c[i,j] * H_i(a) * H_j(b)``
+
+    where the points ``(a,b)`` consist of all pairs of points formed by
+    taking ``a`` from `x` and ``b`` from `y`. The resulting points form a
+    grid with `x` in the first dimension and `y` in the second.
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Hermite series is evaluated at the points in
+        the Cartesian product of `x` and `y`.  If `x` or `y` is a list or
+        tuple, it is first converted to an ndarray, Otherwise it is left
+        unchanged and if it isn't an ndarray it is treated as a scalar. See
+        `hermval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Hermite series at points in the
+        Cartesion product of `x` and `y`.
+
+    See Also
+    --------
+    hermval, hermval2d, hermval3d, hermgrid3d
+
+    """
+    return hermval(y, hermval(x, c))
+
+
+def hermval3d(x, y, z, c):
+    """
+    Evaluate 3D Hermite series at points (x,y,z).
+
+    This function returns the values:
+
+    ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * H_i(x) * H_j(y) * H_k(z)``
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Hermite seres is evaluated at the points
+        ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape.  If
+        any of `x`, `y`, or `z` is a list or tuple, it is first converted
+        to an ndarray. Otherwise it is left unchanged and if it isn't an
+        ndarray it is treated as a scalar.  See `hermval` for explanation of
+        algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Hermite series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    hermval, hermval2d, hermgrid2d, hermgrid3d
+
+    """
+    return hermval(z, hermval2d(x, y, c), False)
+
+
+def hermgrid3d(x, y, z, c):
+    """
+    Evaluate 3D Hermite series on the Cartesian product of x,y,z.
+
+    This function returns the values:
+
+    ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * H_i(a) * H_j(b) * H_k(c)``
+
+    where the points ``(a,b,c)`` consist of all triples formed by taking
+    ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting
+    points form a grid with `x` in the first dimension, `y` in the second,
+    and `z` in the third.
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Hermite seres is evaluated at the points
+        in the cartesian product of `x`, `y`, and `z`
+        ``(x,y,z)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and treated as a scalar.  See
+        `hermval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Hermite series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    hermval, hermval2d, hermgrid2d, hermval3d
+
+    """
+    return hermval(z, hermgrid2d(x, y, c))
+
+
 def hermvander(x, deg) :
     """Vandermonde matrix of given degree.
 
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index be179e47d..17e285e15 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -1,47 +1,51 @@
 """
-Objects for dealing with Hermite series.
+Objects for dealing with Hermite_e series.
 
 This module provides a number of objects (mostly functions) useful for
-dealing with Hermite series, including a `Hermite` class that
+dealing with Hermite_e series, including a `HermiteE` class that
 encapsulates the usual arithmetic operations.  (General information
 on how this module represents and works with such polynomials is in the
 docstring for its "parent" sub-package, `numpy.polynomial`).
 
 Constants
 ---------
-- `hermedomain` -- Hermite series default domain, [-1,1].
-- `hermezero` -- Hermite series that evaluates identically to 0.
-- `hermeone` -- Hermite series that evaluates identically to 1.
-- `hermex` -- Hermite series for the identity map, ``f(x) = x``.
+- `hermedomain` -- Hermite_e series default domain, [-1,1].
+- `hermezero` -- Hermite_e series that evaluates identically to 0.
+- `hermeone` -- Hermite_e series that evaluates identically to 1.
+- `hermex` -- Hermite_e series for the identity map, ``f(x) = x``.
 
 Arithmetic
 ----------
-- `hermemulx` -- multiply a Hermite series in ``P_i(x)`` by ``x``.
-- `hermeadd` -- add two Hermite series.
-- `hermesub` -- subtract one Hermite series from another.
-- `hermemul` -- multiply two Hermite series.
-- `hermediv` -- divide one Hermite series by another.
-- `hermeval` -- evaluate a Hermite series at given points.
+- `hermemulx` -- multiply a Hermite_e series in ``P_i(x)`` by ``x``.
+- `hermeadd` -- add two Hermite_e series.
+- `hermesub` -- subtract one Hermite_e series from another.
+- `hermemul` -- multiply two Hermite_e series.
+- `hermediv` -- divide one Hermite_e series by another.
+- `hermeval` -- evaluate a Hermite_e series at given points.
+- `hermeval2d` -- evaluate a 2D Hermite_e series at given points.
+- `hermeval3d` -- evaluate a 3D Hermite_e series at given points.
+- `hermegrid2d` -- evaluate a 2D Hermite_e series on a Cartesian product.
+- `hermegrid3d` -- evaluate a 3D Hermite_e series on a Cartesian product.
 
 Calculus
 --------
-- `hermeder` -- differentiate a Hermite series.
-- `hermeint` -- integrate a Hermite series.
+- `hermeder` -- differentiate a Hermite_e series.
+- `hermeint` -- integrate a Hermite_e series.
 
 Misc Functions
 --------------
-- `hermefromroots` -- create a Hermite series with specified roots.
-- `hermeroots` -- find the roots of a Hermite series.
-- `hermevander` -- Vandermonde-like matrix for Hermite polynomials.
-- `hermefit` -- least-squares fit returning a Hermite series.
-- `hermetrim` -- trim leading coefficients from a Hermite series.
-- `hermeline` -- Hermite series of given straight line.
-- `herme2poly` -- convert a Hermite series to a polynomial.
-- `poly2herme` -- convert a polynomial to a Hermite series.
+- `hermefromroots` -- create a Hermite_e series with specified roots.
+- `hermeroots` -- find the roots of a Hermite_e series.
+- `hermevander` -- Vandermonde-like matrix for Hermite_e polynomials.
+- `hermefit` -- least-squares fit returning a Hermite_e series.
+- `hermetrim` -- trim leading coefficients from a Hermite_e series.
+- `hermeline` -- Hermite_e series of given straight line.
+- `herme2poly` -- convert a Hermite_e series to a polynomial.
+- `poly2herme` -- convert a polynomial to a Hermite_e series.
 
 Classes
 -------
-- `Hermite` -- A Hermite series class.
+- `HermiteE` -- A Hermite_e series class.
 
 See also
 --------
@@ -57,9 +61,10 @@ import warnings
 from polytemplate import polytemplate
 
 __all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
-        'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermeval',
-        'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots',
-        'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE']
+    'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv', 'hermpow',
+    'legval', 'legval2d', 'legval3d', 'leggrid2d', 'leggrid3d',
+    'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots',
+    'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE']
 
 hermetrim = pu.trimcoef
 
@@ -750,7 +755,7 @@ def hermeint(cs, m=1, k=[], lbnd=0, scl=1):
     >>> from numpy.polynomial.hermite_e import hermeint
     >>> hermeint([1, 2, 3]) # integrate once, value 0 at 0.
     array([ 1.,  1.,  1.,  1.])
-    >>> hermeint([1, 2, 3], m=2) # integrate twice, value & deriv 0 at 0 
+    >>> hermeint([1, 2, 3], m=2) # integrate twice, value & deriv 0 at 0
     array([-0.25      ,  1.        ,  0.5       ,  0.33333333,  0.25      ])
     >>> hermeint([1, 2, 3], k=1) # integrate once, value 1 at 0.
     array([ 2.,  1.,  1.,  1.])
@@ -793,36 +798,58 @@ def hermeint(cs, m=1, k=[], lbnd=0, scl=1):
     return cs
 
 
-def hermeval(x, cs):
-    """Evaluate a Hermite series.
+def hermeval(x, c, tensor=True):
+    """ Evaluate a Hermite_e series.
 
-    If `cs` is of length `n`, this function returns :
+    If `c` is of length ``n + 1``, this function returns the value:
 
-    ``p(x) = cs[0]*P_0(x) + cs[1]*P_1(x) + ... + cs[n-1]*P_{n-1}(x)``
+    ``p(x) = c[0]*He_0(x) + c[1]*He_1(x) + ... + c[n]*He_n(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.
+    If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)``
+    will have the same shape as `x`.  If `x` is a algebra_like object that
+    supports multiplication and addition with itself and the values in `c`,
+    then an object of the same type is returned.
+
+    In the case where c is multidimensional, the shape of the result
+    depends on the value of `tensor`. If tensor is true the shape of the
+    return will be ``c.shape[1:] + x.shape``, where the shape of a scalar
+    is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if
+    `x` is broadcast compatible with that.
+
+    If there are trailing zeros in the coefficients they still take part in
+    the evaluation, so they should be avoided if efficiency is a concern.
 
     Parameters
     ----------
-    x : array_like, ring_like
-        Array of numbers or objects that support multiplication and
-        addition with themselves and with the elements of `cs`.
-    cs : array_like
-        1-d array of Hermite coefficients ordered from low to high.
+    x : array_like, algebra_like
+        If x is a list or tuple, it is converted to an ndarray. Otherwise
+        it is left unchanged and if it isn't an ndarray it is treated as a
+        scalar. In either case, `x` or any element of an ndarray must
+        support addition and multiplication with itself and the elements of
+        `c`.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree n are contained in ``c[n]``. If `c` is multidimesional the
+        remaining indices enumerate multiple Hermite_e series. In the two
+        dimensional case the coefficients may be thought of as stored in
+        the columns of `c`.
+    tensor : boolean, optional
+        If true, the coefficient array shape is extended with ones on the
+        right, one for each dimension of `x`. Scalars are treated as having
+        dimension 0 for this action. The effect is that every column of
+        coefficients in `c` is evaluated for every value in `x`. If False,
+        the `x` are broadcast over the columns of `c` in the usual way.
+        This gives some flexibility in evaluations in the multidimensional
+        case. The default value it ``True``.
 
     Returns
     -------
-    values : ndarray, ring_like
-        If the return is an ndarray then it has the same shape as `x`.
+    values : ndarray, algebra_like
+        The shape of the return value is described above.
 
     See Also
     --------
-    hermefit
-
-    Examples
-    --------
+    hermeval2d, hermegrid2d, hermeval3d, hermegrid3d
 
     Notes
     -----
@@ -839,29 +866,187 @@ def hermeval(x, cs):
            [ 31.,  54.]])
 
     """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    if isinstance(x, tuple) or isinstance(x, list) :
+    c = np.array(c, ndmin=1, copy=0)
+    if c.dtype.char not in 'efdgFDGO':
+        c = c + 0.0
+    if isinstance(x, (tuple, list)):
         x = np.asarray(x)
+    if isinstance(x, np.ndarray) and tensor:
+       c = c.reshape(c.shape + (1,)*x.ndim)
 
-    if len(cs) == 1 :
-        c0 = cs[0]
+    if len(c) == 1 :
+        c0 = c[0]
         c1 = 0
-    elif len(cs) == 2 :
-        c0 = cs[0]
-        c1 = cs[1]
+    elif len(c) == 2 :
+        c0 = c[0]
+        c1 = c[1]
     else :
-        nd = len(cs)
-        c0 = cs[-2]
-        c1 = cs[-1]
-        for i in range(3, len(cs) + 1) :
+        nd = len(c)
+        c0 = c[-2]
+        c1 = c[-1]
+        for i in range(3, len(c) + 1) :
             tmp = c0
             nd =  nd - 1
-            c0 = cs[-i] - c1*(nd - 1)
+            c0 = c[-i] - c1*(nd - 1)
             c1 = tmp + c1*x
     return c0 + c1*x
 
 
+def hermeval2d(x, y, c):
+    """
+    Evaluate 2D Hermite_e series at points (x,y).
+
+    This function returns the values:
+
+    ``p(x,y) = \sum_{i,j} c[i,j] * He_i(x) * He_j(y)``
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Hermite_e seres is evaluated at the points
+        ``(x,y)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray.
+        Otherwise it is left unchanged and if it isn't an ndarray it is
+        treated as a scalar.  See `hermeval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Hermite_e series at points formed
+        from pairs of corresponding values from `x` and `y`.
+
+    See Also
+    --------
+    hermeval, hermegrid2d, hermeval3d, hermegrid3d
+
+    """
+    return hermeval(y, hermeval(x, c), False)
+
+
+def hermegrid2d(x, y, c):
+    """
+    Evaluate 2D Hermite_e series on the Cartesion product of x,y.
+
+    This function returns the values:
+
+    ``p(a,b) = \sum_{i,j} c[i,j] * He_i(a) * He_j(b)``
+
+    where the points ``(a,b)`` consist of all pairs of points formed by
+    taking ``a`` from `x` and ``b`` from `y`. The resulting points form a
+    grid with `x` in the first dimension and `y` in the second.
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Hermite_e series is evaluated at the points in
+        the Cartesian product of `x` and `y`.  If `x` or `y` is a list or
+        tuple, it is first converted to an ndarray, Otherwise it is left
+        unchanged and if it isn't an ndarray it is treated as a scalar. See
+        `hermeval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Hermite_e series at points in the
+        Cartesion product of `x` and `y`.
+
+    See Also
+    --------
+    hermeval, hermeval2d, hermeval3d, hermegrid3d
+
+    """
+    return hermeval(y, hermeval(x, c))
+
+
+def hermeval3d(x, y, z, c):
+    """
+    Evaluate 3D Hermite_e series at points (x,y,z).
+
+    This function returns the values:
+
+    ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * He_i(x) * He_j(y) * He_k(z)``
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Hermite_e seres is evaluated at the points
+        ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape.  If
+        any of `x`, `y`, or `z` is a list or tuple, it is first converted
+        to an ndarray. Otherwise it is left unchanged and if it isn't an
+        ndarray it is treated as a scalar.  See `hermeval` for explanation of
+        algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Hermite_e series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    hermeval, hermeval2d, hermegrid2d, hermegrid3d
+
+    """
+    return hermeval(z, hermeval2d(x, y, c), False)
+
+
+def hermegrid3d(x, y, z, c):
+    """
+    Evaluate 3D Hermite_e series on the Cartesian product of x,y,z.
+
+    This function returns the values:
+
+    ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * He_i(a) * He_j(b) * He_k(c)``
+
+    where the points ``(a,b,c)`` consist of all triples formed by taking
+    ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting
+    points form a grid with `x` in the first dimension, `y` in the second,
+    and `z` in the third.
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Hermite_e seres is evaluated at the points
+        in the cartesian product of `x`, `y`, and `z`
+        ``(x,y,z)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and treated as a scalar.  See
+        `hermeval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Hermite_e series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    hermeval, hermeval2d, hermegrid2d, hermeval3d
+
+    """
+    return hermeval(z, hermegrid2d(x, y, c))
+
+
 def hermevander(x, deg) :
     """Vandermonde matrix of given degree.
 
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index e33707895..dbbb2c401 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -22,6 +22,10 @@ Arithmetic
 - `lagmul` -- multiply two Laguerre series.
 - `lagdiv` -- divide one Laguerre series by another.
 - `lagval` -- evaluate a Laguerre series at given points.
+- `lagval2d` -- evaluate a 2D Laguerre series at given points.
+- `lagval3d` -- evaluate a 3D Laguerre series at given points.
+- `laggrid2d` -- evaluate a 2D Laguerre series on a Cartesian product.
+- `laggrid3d` -- evaluate a 3D Laguerre series on a Cartesian product.
 
 Calculus
 --------
@@ -50,17 +54,18 @@ See also
 """
 from __future__ import division
 
-__all__ = ['lagzero', 'lagone', 'lagx', 'lagdomain', 'lagline',
-        'lagadd', 'lagsub', 'lagmulx', 'lagmul', 'lagdiv', 'lagval',
-        'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots',
-        'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre']
-
 import numpy as np
 import numpy.linalg as la
 import polyutils as pu
 import warnings
 from polytemplate import polytemplate
 
+__all__ = ['lagzero', 'lagone', 'lagx', 'lagdomain', 'lagline',
+    'lagadd', 'lagsub', 'lagmulx', 'lagmul', 'lagdiv', 'lagpow',
+    'lagval', 'lagval2d', 'lagval3d', 'laggrid2d', 'laggrid3d',
+    'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots',
+    'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre']
+
 lagtrim = pu.trimcoef
 
 
@@ -794,36 +799,59 @@ def lagint(cs, m=1, k=[], lbnd=0, scl=1):
     return cs
 
 
-def lagval(x, cs):
-    """Evaluate a Laguerre series.
+def lagval(x, c, tensor=True):
+    """
+    Evaluate a Laguerre series.
+
+    If `c` is of length ``n + 1``, this function returns the value:
 
-    If `cs` is of length `n`, this function returns :
+    ``p(x) = c[0]*L_0(x) + c[1]*L_1(x) + ... + c[n]*L_n(x)``
 
-    ``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 and `c` is 1 dimensional, then ``p(x)``
+    will have the same shape as `x`.  If `x` is a algebra_like object that
+    supports multiplication and addition with itself and the values in `c`,
+    then an object of the same type is returned.
 
-    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.
+    In the case where c is multidimensional, the shape of the result
+    depends on the value of `tensor`. If tensor is true the shape of the
+    return will be ``c.shape[1:] + x.shape``, where the shape of a scalar
+    is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if
+    `x` is broadcast compatible with that.
+
+    If there are trailing zeros in the coefficients they still take part in
+    the evaluation, so they should be avoided if efficiency is a concern.
 
     Parameters
     ----------
-    x : array_like, ring_like
-        Array of numbers or objects that support multiplication and
-        addition with themselves and with the elements of `cs`.
-    cs : array_like
-        1-d array of Laguerre coefficients ordered from low to high.
+    x : array_like, algebra_like
+        If x is a list or tuple, it is converted to an ndarray. Otherwise
+        it is left unchanged and if it isn't an ndarray it is treated as a
+        scalar. In either case, `x` or any element of an ndarray must
+        support addition and multiplication with itself and the elements of
+        `c`.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree n are contained in ``c[n]``. If `c` is multidimesional the
+        remaining indices enumerate multiple Laguerre series. In the two
+        dimensional case the coefficients may be thought of as stored in
+        the columns of `c`.
+    tensor : boolean, optional
+        If true, the coefficient array shape is extended with ones on the
+        right, one for each dimension of `x`. Scalars are treated as having
+        dimension 0 for this action. The effect is that every column of
+        coefficients in `c` is evaluated for every value in `x`. If False,
+        the `x` are broadcast over the columns of `c` in the usual way.
+        This gives some flexibility in evaluations in the multidimensional
+        case. The default value it ``True``.
 
     Returns
     -------
-    values : ndarray, ring_like
-        If the return is an ndarray then it has the same shape as `x`.
+    values : ndarray, algebra_like
+        The shape of the return value is described above.
 
     See Also
     --------
-    lagfit
-
-    Examples
-    --------
+    lagval2d, laggrid2d, lagval3d, laggrid3d
 
     Notes
     -----
@@ -840,29 +868,188 @@ def lagval(x, cs):
            [-4.5, -2. ]])
 
     """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    if isinstance(x, tuple) or isinstance(x, list) :
+    c = np.array(c, ndmin=1, copy=0)
+    if c.dtype.char not in 'efdgFDGO':
+        c = c + 0.0
+    if isinstance(x, (tuple, list)):
         x = np.asarray(x)
+    if isinstance(x, np.ndarray) and tensor:
+       c = c.reshape(c.shape + (1,)*x.ndim)
+
 
-    if len(cs) == 1 :
-        c0 = cs[0]
+    if len(c) == 1 :
+        c0 = c[0]
         c1 = 0
-    elif len(cs) == 2 :
-        c0 = cs[0]
-        c1 = cs[1]
+    elif len(c) == 2 :
+        c0 = c[0]
+        c1 = c[1]
     else :
-        nd = len(cs)
-        c0 = cs[-2]
-        c1 = cs[-1]
-        for i in range(3, len(cs) + 1) :
+        nd = len(c)
+        c0 = c[-2]
+        c1 = c[-1]
+        for i in range(3, len(c) + 1) :
             tmp = c0
             nd =  nd - 1
-            c0 = cs[-i] - (c1*(nd - 1))/nd
+            c0 = c[-i] - (c1*(nd - 1))/nd
             c1 = tmp + (c1*((2*nd - 1) - x))/nd
     return c0 + c1*(1 - x)
 
 
+def lagval2d(x, y, c):
+    """
+    Evaluate 2D Laguerre series at points (x,y).
+
+    This function returns the values:
+
+    ``p(x,y) = \sum_{i,j} c[i,j] * L_i(x) * L_j(y)``
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Laguerre seres is evaluated at the points
+        ``(x,y)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray.
+        Otherwise it is left unchanged and if it isn't an ndarray it is
+        treated as a scalar.  See `lagval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Laguerre series at points formed
+        from pairs of corresponding values from `x` and `y`.
+
+    See Also
+    --------
+    lagval, laggrid2d, lagval3d, laggrid3d
+
+    """
+    return lagval(y, lagval(x, c), False)
+
+
+def laggrid2d(x, y, c):
+    """
+    Evaluate 2D Laguerre series on the Cartesion product of x,y.
+
+    This function returns the values:
+
+    ``p(a,b) = \sum_{i,j} c[i,j] * L_i(a) * L_j(b)``
+
+    where the points ``(a,b)`` consist of all pairs of points formed by
+    taking ``a`` from `x` and ``b`` from `y`. The resulting points form a
+    grid with `x` in the first dimension and `y` in the second.
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Laguerre series is evaluated at the points in
+        the Cartesian product of `x` and `y`.  If `x` or `y` is a list or
+        tuple, it is first converted to an ndarray, Otherwise it is left
+        unchanged and if it isn't an ndarray it is treated as a scalar. See
+        `lagval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Laguerre series at points in the
+        Cartesion product of `x` and `y`.
+
+    See Also
+    --------
+    lagval, lagval2d, lagval3d, laggrid3d
+
+    """
+    return lagval(y, lagval(x, c))
+
+
+def lagval3d(x, y, z, c):
+    """
+    Evaluate 3D Laguerre series at points (x,y,z).
+
+    This function returns the values:
+
+    ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * L_i(x) * L_j(y) * L_k(z)``
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Laguerre seres is evaluated at the points
+        ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape.  If
+        any of `x`, `y`, or `z` is a list or tuple, it is first converted
+        to an ndarray. Otherwise it is left unchanged and if it isn't an
+        ndarray it is treated as a scalar.  See `lagval` for explanation of
+        algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Laguerre series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    lagval, lagval2d, laggrid2d, laggrid3d
+
+    """
+    return lagval(z, lagval2d(x, y, c), False)
+
+
+def laggrid3d(x, y, z, c):
+    """
+    Evaluate 3D Laguerre series on the Cartesian product of x,y,z.
+
+    This function returns the values:
+
+    ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * L_i(a) * L_j(b) * L_k(c)``
+
+    where the points ``(a,b,c)`` consist of all triples formed by taking
+    ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting
+    points form a grid with `x` in the first dimension, `y` in the second,
+    and `z` in the third.
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Laguerre seres is evaluated at the points
+        in the cartesian product of `x`, `y`, and `z`
+        ``(x,y,z)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and treated as a scalar.  See
+        `lagval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Laguerre series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    lagval, lagval2d, laggrid2d, lagval3d
+
+    """
+    return lagval(z, laggrid2d(x, y, c))
+
+
 def lagvander(x, deg) :
     """Vandermonde matrix of given degree.
 
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index 4a2082e3c..5a72217d0 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -23,6 +23,10 @@ Arithmetic
 - `legdiv` -- divide one Legendre series by another.
 - `legpow` -- raise a Legendre series to an positive integer power
 - `legval` -- evaluate a Legendre series at given points.
+- `legval2d` -- evaluate a 2D Legendre series at given points.
+- `legval3d` -- evaluate a 3D Legendre series at given points.
+- `leggrid2d` -- evaluate a 2D Legendre series on a Cartesian product.
+- `leggrid3d` -- evaluate a 3D Legendre series on a Cartesian product.
 
 Calculus
 --------
@@ -51,18 +55,18 @@ See also
 """
 from __future__ import division
 
-__all__ = ['legzero', 'legone', 'legx', 'legdomain', 'legline',
-        'legadd', 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow',
-        '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
 
+__all__ = ['legzero', 'legone', 'legx', 'legdomain', 'legline',
+    'legadd', 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow',
+    'legval', 'legval2d', 'legval3d', 'leggrid2d', 'leggrid3d',
+    'legder', 'legint', 'leg2poly', 'poly2leg', 'legfromroots',
+    'legvander', 'legfit', 'legtrim', 'legroots', 'Legendre']
+
 legtrim = pu.trimcoef
 
 
@@ -814,33 +818,58 @@ def legint(cs, m=1, k=[], lbnd=0, scl=1):
     return cs
 
 
-def legval(x, cs):
-    """Evaluate a Legendre series.
+def legval(x, c, tensor=True):
+    """ Evaluate a Legendre series.
+
+    If `c` is of length ``n + 1``, this function returns the value:
 
-    If `cs` is of length `n`, this function returns :
+    ``p(x) = c[0]*L_0(x) + c[1]*L_1(x) + ... + c[n]*L_n(x)``
 
-    ``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 and `c` is 1 dimensional, then ``p(x)``
+    will have the same shape as `x`.  If `x` is a algebra_like object that
+    supports multiplication and addition with itself and the values in `c`,
+    then an object of the same type is returned.
 
-    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.
+    In the case where c is multidimensional, the shape of the result
+    depends on the value of `tensor`. If tensor is true the shape of the
+    return will be ``c.shape[1:] + x.shape``, where the shape of a scalar
+    is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if
+    `x` is broadcast compatible with that.
+
+    If there are trailing zeros in the coefficients they still take part in
+    the evaluation, so they should be avoided if efficiency is a concern.
 
     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 Legendre coefficients ordered from low to high.
+    x : array_like, algebra_like
+        If x is a list or tuple, it is converted to an ndarray. Otherwise
+        it is left unchanged and if it isn't an ndarray it is treated as a
+        scalar. In either case, `x` or any element of an ndarray must
+        support addition and multiplication with itself and the elements of
+        `c`.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree n are contained in ``c[n]``. If `c` is multidimesional the
+        remaining indices enumerate multiple Legendre series. In the two
+        dimensional case the coefficients may be thought of as stored in
+        the columns of `c`.
+    tensor : boolean, optional
+        If true, the coefficient array shape is extended with ones on the
+        right, one for each dimension of `x`. Scalars are treated as having
+        dimension 0 for this action. The effect is that every column of
+        coefficients in `c` is evaluated for every value in `x`. If False,
+        the `x` are broadcast over the columns of `c` in the usual way.
+        This gives some flexibility in evaluations in the multidimensional
+        case. The default value it ``True``.
 
     Returns
     -------
-    values : ndarray, ring_like
-        If the return is an ndarray then it has the same shape as `x`.
+    values : ndarray, algebra_like
+        The shape of the return value is described above.
 
     See Also
     --------
-    legfit
+    legval2d, leggrid2d, legval3d, leggrid3d
 
     Notes
     -----
@@ -850,29 +879,187 @@ def legval(x, cs):
     --------
 
     """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    if isinstance(x, tuple) or isinstance(x, list) :
+    c = np.array(c, ndmin=1, copy=0)
+    if c.dtype.char not in 'efdgFDGO':
+        c = c + 0.0
+    if isinstance(x, (tuple, list)):
         x = np.asarray(x)
+    if isinstance(x, np.ndarray) and tensor:
+       c = c.reshape(c.shape + (1,)*x.ndim)
 
-    if len(cs) == 1 :
-        c0 = cs[0]
+    if len(c) == 1 :
+        c0 = c[0]
         c1 = 0
-    elif len(cs) == 2 :
-        c0 = cs[0]
-        c1 = cs[1]
+    elif len(c) == 2 :
+        c0 = c[0]
+        c1 = c[1]
     else :
-        nd = len(cs)
-        c0 = cs[-2]
-        c1 = cs[-1]
-        for i in range(3, len(cs) + 1) :
+        nd = len(c)
+        c0 = c[-2]
+        c1 = c[-1]
+        for i in range(3, len(c) + 1) :
             tmp = c0
             nd =  nd - 1
-            c0 = cs[-i] - (c1*(nd - 1))/nd
+            c0 = c[-i] - (c1*(nd - 1))/nd
             c1 = tmp + (c1*x*(2*nd - 1))/nd
     return c0 + c1*x
 
 
+def legval2d(x, y, c):
+    """
+    Evaluate 2D Legendre series at points (x,y).
+
+    This function returns the values:
+
+    ``p(x,y) = \sum_{i,j} c[i,j] * L_i(x) * L_j(y)``
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Legendre seres is evaluated at the points
+        ``(x,y)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray.
+        Otherwise it is left unchanged and if it isn't an ndarray it is
+        treated as a scalar.  See `legval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Legendre series at points formed
+        from pairs of corresponding values from `x` and `y`.
+
+    See Also
+    --------
+    legval, leggrid2d, legval3d, leggrid3d
+
+    """
+    return legval(y, legval(x, c), False)
+
+
+def leggrid2d(x, y, c):
+    """
+    Evaluate 2D Legendre series on the Cartesion product of x,y.
+
+    This function returns the values:
+
+    ``p(a,b) = \sum_{i,j} c[i,j] * L_i(a) * L_j(b)``
+
+    where the points ``(a,b)`` consist of all pairs of points formed by
+    taking ``a`` from `x` and ``b`` from `y`. The resulting points form a
+    grid with `x` in the first dimension and `y` in the second.
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional Legendre series is evaluated at the points in
+        the Cartesian product of `x` and `y`.  If `x` or `y` is a list or
+        tuple, it is first converted to an ndarray, Otherwise it is left
+        unchanged and if it isn't an ndarray it is treated as a scalar. See
+        `legval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional Legendre series at points in the
+        Cartesion product of `x` and `y`.
+
+    See Also
+    --------
+    legval, legval2d, legval3d, leggrid3d
+
+    """
+    return legval(y, legval(x, c))
+
+
+def legval3d(x, y, z, c):
+    """
+    Evaluate 3D Legendre series at points (x,y,z).
+
+    This function returns the values:
+
+    ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * L_i(x) * L_j(y) * L_k(z)``
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Legendre seres is evaluated at the points
+        ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape.  If
+        any of `x`, `y`, or `z` is a list or tuple, it is first converted
+        to an ndarray. Otherwise it is left unchanged and if it isn't an
+        ndarray it is treated as a scalar.  See `legval` for explanation of
+        algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Legendre series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    legval, legval2d, leggrid2d, leggrid3d
+
+    """
+    return legval(z, legval2d(x, y, c), False)
+
+
+def leggrid3d(x, y, z, c):
+    """
+    Evaluate 3D Legendre series on the Cartesian product of x,y,z.
+
+    This function returns the values:
+
+    ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * L_i(a) * L_j(b) * L_k(c)``
+
+    where the points ``(a,b,c)`` consist of all triples formed by taking
+    ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting
+    points form a grid with `x` in the first dimension, `y` in the second,
+    and `z` in the third.
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional Legendre seres is evaluated at the points
+        in the cartesian product of `x`, `y`, and `z`
+        ``(x,y,z)``, where `x` and `y` must have the same shape.  If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and treated as a scalar.  See
+        `legval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the three dimensional Legendre series at points formed
+        from triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    legval, legval2d, leggrid2d, legval3d
+
+    """
+    return legval(z, leggrid2d(x, y, c))
+
+
 def legvander(x, deg) :
     """Vandermonde matrix of given degree.
 
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index a0716b9a5..a0c9f903a 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -22,6 +22,10 @@ Arithmetic
 - `polydiv` -- divide one polynomial by another.
 - `polypow` -- raise a polynomial to an positive integer power
 - `polyval` -- evaluate a polynomial at given points.
+- `polyval2d` -- evaluate a 2D polynomial at given points.
+- `polyval3d` -- evaluate a 3D polynomial at given points.
+- `polygrid2d` -- evaluate a 2D polynomial on a Cartesian product.
+- `polygrid3d` -- evaluate a 3D polynomial on a Cartesian product.
 
 Calculus
 --------
@@ -50,8 +54,9 @@ from __future__ import division
 
 __all__ = ['polyzero', 'polyone', 'polyx', 'polydomain', 'polyline',
     'polyadd', 'polysub', 'polymulx', 'polymul', 'polydiv', 'polypow',
-    'polyval', 'polyder', 'polyint', 'polyfromroots', 'polyvander',
-    'polyfit', 'polytrim', 'polyroots', 'Polynomial']
+    'polyval', 'polyval2d', 'polyval3d', 'polygrid2d', 'polygrid3d',
+    'polyder', 'polyint', 'polyfromroots', 'polyvander', 'polyfit',
+    'polytrim', 'polyroots', 'Polynomial']
 
 import numpy as np
 import numpy.linalg as la
@@ -615,52 +620,255 @@ def polyint(cs, m=1, k=[], lbnd=0, scl=1):
     return cs
 
 
-def polyval(x, cs):
+def polyval(x, c, tensor=True):
     """
     Evaluate a polynomial.
 
-    If `cs` is of length `n`, this function returns :
+    If `c` is of length ``n + 1``, this function returns the value:
 
-    ``p(x) = cs[0] + cs[1]*x + ... + cs[n-1]*x**(n-1)``
+    ``p(x) = c[0] + c[1]*x + ... + c[n]*x**(n)``
 
-    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.
+    If `x` is a sequence or array and `c` is 1 dimensional, then ``p(x)``
+    will have the same shape as `x`.  If `x` is a algebra_like object that
+    supports multiplication and addition with itself and the values in `c`,
+    then an object of the same type is returned.
+
+    In the case where c is multidimensional, the shape of the result
+    depends on the value of `tensor`. If tensor is true the shape of the
+    return will be ``c.shape[1:] + x.shape``, where the shape of a scalar
+    is the empty tuple. If tensor is false the shape is ``c.shape[1:]`` if
+    `x` is broadcast compatible with that.
+
+    If there are trailing zeros in the coefficients they still take part in
+    the evaluation, so they should be avoided if efficiency is a concern.
 
     Parameters
     ----------
-    x : array_like, ring_like
+    x : array_like, algebra_like
         If x is a list or tuple, it is converted to an ndarray. Otherwise
-        it must support addition and multiplication with itself and the
-        elements of `cs`.
-    cs : array_like
-        1-d array of Chebyshev coefficients ordered from low to high.
+        it is left unchanged and if it isn't an ndarray it is treated as a
+        scalar. In either case, `x` or any element of an ndarray must
+        support addition and multiplication with itself and the elements of
+        `c`.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree n are contained in ``c[n]``. If `c` is multidimesional the
+        remaining indices enumerate multiple polynomials. In the two
+        dimensional case the coefficients may be thought of as stored in
+        the columns of `c`.
+    tensor : boolean, optional
+        If true, the coefficient array shape is extended with ones on the
+        right, one for each dimension of `x`. Scalars are treated as having
+        dimension 0 for this action. The effect is that every column of
+        coefficients in `c` is evaluated for every value in `x`. If False,
+        the `x` are broadcast over the columns of `c` in the usual way.
+        This gives some flexibility in evaluations in the multidimensional
+        case. The default value it ``True``.
 
     Returns
     -------
-    values : ndarray
-        The return array has the same shape as `x`.
+    values : ndarray, algebra_like
+        The shape of the return value is described above.
 
     See Also
     --------
-    polyfit
+    polyval2d, polygrid2d, polyval3d, polygrid3d
 
     Notes
     -----
     The evaluation uses Horner's method.
 
+    Examples
+    --------
+    >>> from numpy.polynomial.polynomial import polyval
+    >>> polyval(1, [1,2,3])
+    6.0
+    >>> a = np.arange(4).reshape(2,2)
+    >>> a
+    array([[[  0,   1],
+            [  2,   3]],
+    >>> polyval(a, [1,2,3])
+    array([[  1.,   6.],
+           [ 17.,  34.]])
+    >>> c = np.arange(4).reshape(2,2)
+    >>> c
+    array([[[  0,   1],
+            [  2,   3]],
+    >>> polyval([1,2], c, tensor=True)
+    array([[ 2.,  4.],
+           [ 4.,  7.]])
+    >>> polyval([1,2], c, tensor=False)
+    array([ 2.,  7.])
+
     """
-    # cs is a trimmed copy
-    [cs] = pu.as_series([cs])
-    if isinstance(x, tuple) or isinstance(x, list) :
+    c = np.array(c, ndmin=1, copy=0)
+    if c.dtype.char not in 'efdgFDGO':
+        c = c + 0.0
+    if isinstance(x, (tuple, list)):
         x = np.asarray(x)
+    if isinstance(x, np.ndarray) and tensor:
+       c = c.reshape(c.shape + (1,)*x.ndim)
 
-    c0 = cs[-1] + x*0
-    for i in range(2, len(cs) + 1) :
-        c0 = cs[-i] + c0*x
+    c0 = c[-1] + x*0
+    for i in range(2, len(c) + 1) :
+        c0 = c[-i] + c0*x
     return c0
 
 
+def polyval2d(x, y, c):
+    """
+    Evaluate 2D polynomials at points (x,y).
+
+    This function returns the values:
+
+    ``p(x,y) = \sum_{i,j} c[i,j] * x^i * y^j``
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional polynomial is evaluated at the points
+        ``(x,y)``, where `x` and `y` must have the same shape. If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and if it isn't an ndarray it is
+        treated as a scalar. See `polyval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional polynomial at points formed with
+        pairs of corresponding values from `x` and `y`.
+
+    See Also
+    --------
+    polyval, polygrid2d, polyval3d, polygrid3d
+
+    """
+    return polyval(y, polyval(x, c), False)
+
+
+def polygrid2d(x, y, c):
+    """
+    Evaluate 2D polynomials on the Cartesion product of x,y.
+
+    This function returns the values:
+
+    ``p(a,b) = \sum_{i,j} c[i,j] * a^i * b^j``
+
+    where the points ``(a,b)`` consist of all pairs of points formed by
+    taking ``a`` from `x` and ``b`` from `y`. The resulting points form a
+    grid with `x` in the first dimension and `y` in the second.
+
+    Parameters
+    ----------
+    x,y : array_like, algebra_like
+        The two dimensional polynomial is evaluated at the points in the
+        Cartesian product of `x` and `y`.  If `x` or `y` is a list or
+        tuple, it is first converted to an ndarray, otherwise it is left
+        unchanged and if it isn't an ndarray it is treated as a scalar. See
+        `polyval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j]``. If `c` has dimension
+        greater than 2 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the two dimensional polynomial at points in the Cartesion
+        product of `x` and `y`.
+
+    See Also
+    --------
+    polyval, polyval2d, polyval3d, polygrid3d
+
+    """
+    return polyval(y, polyval(x, c))
+
+
+def polyval3d(x, y, z, c):
+    """
+    Evaluate 3D polynomials at points (x,y,z).
+
+    This function returns the values:
+
+    ``p(x,y,z) = \sum_{i,j,k} c[i,j,k] * x^i * y^j * z^k``
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional polynomial is evaluated at the points
+        ``(x,y,z)``, where `x`, `y`, and `z` must have the same shape.  If
+        any of `x`, `y`, or `z` is a list or tuple, it is first converted
+        to an ndarray, otherwise it is left unchanged and if it isn't an
+        ndarray it is  treated as a scalar. See `polyval` for an
+        explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i, j are contained in ``c[i,j,k]``. If `c` has dimension
+        greater than 3 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the multidimension polynomial on points formed with
+        triples of corresponding values from `x`, `y`, and `z`.
+
+    See Also
+    --------
+    polyval, polyval2d, polygrid2d, polygrid3d
+
+    """
+    return polyval(z, polyval2d(x, y, c), False)
+
+
+def polygrid3d(x, y, z, c):
+    """
+    Evaluate 3D polynomials on the Cartesion product of x,y,z.
+
+    This function returns the values:
+
+    ``p(a,b,c) = \sum_{i,j,k} c[i,j,k] * a^i * b^j * c^k``
+
+    where the points ``(a,b,c)`` consist of all triples formed by taking
+    ``a`` from `x`, ``b`` from `y`, and ``c`` from `z`. The resulting
+    points form a grid with `x` in the first dimension, `y` in the second,
+    and `z` in the third.
+
+    Parameters
+    ----------
+    x,y,z : array_like, algebra_like
+        The three dimensional polynomial is evaluated at the points
+        ``(x,y,z)``, where `x` and `y` must have the same shape. If `x` or
+        `y` is a list or tuple, it is first converted to an ndarray,
+        otherwise it is left unchanged and treated as a scalar. See
+        `polyval` for explanation of algebra_like.
+    c : array_like
+        Array of coefficients ordered so that the coefficients for terms of
+        degree i,j are contained in ``c[i,j,k]``. If `c` has dimension
+        greater than 3 the remaining indices enumerate multiple sets of
+        coefficients.
+
+    Returns
+    -------
+    values : ndarray, algebra_like
+        The values of the multidimensional polynomial on the cartesion
+        product of `x`, `y`, and `z`.
+
+    See Also
+    --------
+    polyval, polyval2d, polygrid2d, polyval3d
+
+    """
+    return polyval(z, polygrid2d(x, y, c))
+
+
 def polyvander(x, deg) :
     """Vandermonde matrix of given degree.