ENH: Add functions for Gauss quadrature and associated weight functions.

The new functions for Gauss quadrature  are of the form xxxgauss, where xxx
is any of cheb, leg, lag, herm, herme. They return the Gauss points and
weights for Gauss quadrature of the various orthogonal polynomial types
given the degree. They are tested to work up to degree 100.

The new functions for the weight are of the form xxxweight, where xxx is
any of cheb, leg, lag, herm, herme. They return the value of the weight
function for the various orthogonal polynomial types given and array of
points.

5 files changed, 433 insertions, 5 deletions

diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index c80123bd3..dbadaa67b 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -42,6 +42,9 @@ Misc Functions
 - `chebvander` -- Vandermonde-like matrix for Chebyshev polynomials.
 - `chebvander2d` -- Vandermonde-like matrix for 2D power series.
 - `chebvander3d` -- Vandermonde-like matrix for 3D power series.
+- `chebgauss` -- Gauss-Chebyshev quadrature, points and weights.
+- `chebweight` -- Chebyshev weight function.
+- `chebcompanion` -- symmetrized companion matrix in Chebyshev form.
 - `chebfit` -- least-squares fit returning a Chebyshev series.
 - `chebpts1` -- Chebyshev points of the first kind.
 - `chebpts2` -- Chebyshev points of the second kind.
@@ -95,7 +98,8 @@ __all__ = ['chebzero', 'chebone', 'chebx', 'chebdomain', 'chebline',
     'chebval', 'chebder', 'chebint', 'cheb2poly', 'poly2cheb',
     'chebfromroots', 'chebvander', 'chebfit', 'chebtrim', 'chebroots',
     'chebpts1', 'chebpts2', 'Chebyshev', 'chebval2d', 'chebval3d',
-    'chebgrid2d', 'chebgrid3d', 'chebvander2d','chebvander3d']
+    'chebgrid2d', 'chebgrid3d', 'chebvander2d','chebvander3d',
+    'chebcompanion', 'chebgauss', 'chebweight']
 
 chebtrim = pu.trimcoef
 
@@ -1724,6 +1728,72 @@ def chebroots(cs):
     return r
 
 
+def chebgauss(deg):
+    """Gauss Chebyshev quadrature.
+
+    Computes the sample points and weights for Gauss-Chebyshev quadrature.
+    These sample points and weights will correctly integrate polynomials of
+    degree ``2*deg - 1`` or less over the interval ``[-1, 1]`` with the
+    weight function ``f(x) = 1/sqrt(1 - x**2)``.
+
+    Parameters
+    ----------
+    deg : int
+        Number of sample points and weights. It must be >= 1.
+
+    Returns
+    -------
+    x : ndarray
+        1-D ndarray containing the sample points.
+    y : ndarray
+        1-D ndarray containing the weights.
+
+    Notes
+    -----
+    The results have only been tested up to degree 100. Higher degrees may
+    be problematic. There are closed form solutions for the sample points
+    and weights. If ``n = deg``, then
+
+          ``x_i = cos(pi*(2*i - 1)/(2*n))``
+          ``w_i = pi/n``
+
+    where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th
+    root of ``L_n``, and then scaling the results to get the right value
+    when integrating 1.
+
+    """
+    ideg = int(deg)
+    if ideg != deg or ideg < 1:
+        raise ValueError("deg must be a non-negative integer")
+
+    x = np.cos(np.pi * np.arange(1, 2*ideg, 2) / (2.0*ideg))
+    w = np.ones(ideg)*(np.pi/ideg)
+
+    return x, w
+
+
+def chebweight(x):
+    """Weight function of the Chebyshev polynomials.
+
+    The weight function for which the Chebyshev polynomials are orthogonal.
+    In this case the weight function is ``1/(1 - x**2)``. Note that the
+    Chebyshev polynomials are not normalized.
+
+    Parameters
+    ----------
+    x : array_like
+       Values at which the weight function will be computed.
+
+    Returns
+    -------
+    w : ndarray
+       The weight function at `x`.
+
+    """
+    w = 1./(np.sqrt(1. + x) * np.sqrt(1. - x))
+    return w
+
+
 def chebpts1(npts):
     """Chebyshev points of the first kind.
 
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 5babfb7b1..58cc655d7 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -39,6 +39,9 @@ Misc Functions
 - `hermvander` -- Vandermonde-like matrix for Hermite polynomials.
 - `hermvander2d` -- Vandermonde-like matrix for 2D power series.
 - `hermvander3d` -- Vandermonde-like matrix for 3D power series.
+- `hermgauss` -- Gauss-Hermite quadrature, points and weights.
+- `hermweight` -- Hermite weight function.
+- `hermcompanion` -- symmetrized companion matrix in Hermite form.
 - `hermfit` -- least-squares fit returning a Hermite series.
 - `hermtrim` -- trim leading coefficients from a Hermite series.
 - `hermline` -- Hermite series of given straight line.
@@ -67,7 +70,8 @@ __all__ = ['hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline',
     'hermval', 'hermder', 'hermint', 'herm2poly', 'poly2herm',
     'hermfromroots', 'hermvander', 'hermfit', 'hermtrim', 'hermroots',
     'Hermite', 'hermval2d', 'hermval3d', 'hermgrid2d', 'hermgrid3d',
-    'hermvander2d', 'hermvander3d']
+    'hermvander2d', 'hermvander3d', 'hermcompanion', 'hermgauss',
+    'hermweight']
 
 hermtrim = pu.trimcoef
 
@@ -1507,6 +1511,93 @@ def hermroots(cs):
     return r
 
 
+def hermgauss(deg):
+    """Gauss Hermite quadrature.
+
+    Computes the sample points and weights for Gauss-Hermite quadrature.
+    These sample points and weights will correctly integrate polynomials of
+    degree ``2*deg - 1`` or less over the interval ``[-inf, inf]`` with the
+    weight function ``f(x) = exp(-x**2)``.
+
+    Parameters
+    ----------
+    deg : int
+        Number of sample points and weights. It must be >= 1.
+
+    Returns
+    -------
+    x : ndarray
+        1-D ndarray containing the sample points.
+    y : ndarray
+        1-D ndarray containing the weights.
+
+    Notes
+    -----
+    The results have only been tested up to degree 100. Higher degrees may
+    be problematic. The weights are determined by using the fact that
+
+          w = c / (H'_n(x_k) * H_{n-1}(x_k))
+
+    where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th
+    root of ``H_n``, and then scaling the results to get the right value
+    when integrating 1.
+
+    """
+    ideg = int(deg)
+    if ideg != deg or ideg < 1:
+        raise ValueError("deg must be a non-negative integer")
+
+    # first approximation of roots. We use the fact that the companion
+    # matrix is symmetric in this case in order to obtain better zeros.
+    c = np.array([0]*deg + [1])
+    m = hermcompanion(c)
+    x = la.eigvals(m)
+    x.sort()
+
+    # improve roots by one application of Newton
+    dy = hermval(x, c)
+    df = hermval(x, hermder(c))
+    x -= dy/df
+
+    # compute the weights. We scale the factor to avoid possible numerical
+    # overflow.
+    fm = hermval(x, c[1:])
+    fm /= np.abs(fm).max()
+    df /= np.abs(df).max()
+    w = 1/(fm * df)
+
+    # for Hermite we can also symmetrize
+    w = (w + w[::-1])/2
+    x = (x - x[::-1])/2
+
+    # scale w to get the right value
+    w *= np.sqrt(np.pi) / w.sum()
+
+    return x, w
+
+
+def hermweight(x):
+    """Weight function of the Hermite polynomials.
+
+    The weight function for which the Hermite polynomials are orthogonal.
+    In this case the weight function is ``exp(-x**2)``. Note that the
+    Hermite polynomials are not normalized.
+
+    Parameters
+    ----------
+    x : array_like
+       Values at which the weight function will be computed.
+
+    Returns
+    -------
+    w : ndarray
+       The weight function at `x`.
+
+    """
+    w = np.exp(-x**2)
+    return w
+
+
 #
 # Hermite series class
 #
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index 97ef366bc..c8d641734 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -39,6 +39,9 @@ Misc Functions
 - `hermevander` -- Vandermonde-like matrix for Hermite_e polynomials.
 - `hermevander2d` -- Vandermonde-like matrix for 2D power series.
 - `hermevander3d` -- Vandermonde-like matrix for 3D power series.
+- `hermegauss` -- Gauss-Hermite_e quadrature, points and weights.
+- `hermeweight` -- Hermite_e weight function.
+- `hermecompanion` -- symmetrized companion matrix in Hermite_e form.
 - `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.
@@ -68,7 +71,7 @@ __all__ = ['hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
     'hermeder', 'hermeint', 'herme2poly', 'poly2herme', 'hermefromroots',
     'hermevander', 'hermefit', 'hermetrim', 'hermeroots', 'HermiteE',
     'hermeval2d', 'hermeval3d', 'hermegrid2d', 'hermegrid3d', 'hermevander2d',
-    'hermevander3d']
+    'hermevander3d', 'hermecompanion', 'hermegauss', 'hermeweight']
 
 hermetrim = pu.trimcoef
 
@@ -1503,6 +1506,93 @@ def hermeroots(cs):
     return r
 
 
+def hermegauss(deg):
+    """Gauss Hermite_e quadrature.
+
+    Computes the sample points and weights for Gauss-Hermite_e quadrature.
+    These sample points and weights will correctly integrate polynomials of
+    degree ``2*deg - 1`` or less over the interval ``[-inf, inf]`` with the
+    weight function ``f(x) = exp(-.5*x**2)``.
+
+    Parameters
+    ----------
+    deg : int
+        Number of sample points and weights. It must be >= 1.
+
+    Returns
+    -------
+    x : ndarray
+        1-D ndarray containing the sample points.
+    y : ndarray
+        1-D ndarray containing the weights.
+
+    Notes
+    -----
+    The results have only been tested up to degree 100. Higher degrees may
+    be problematic. The weights are determined by using the fact that
+
+          w = c / (He'_n(x_k) * He_{n-1}(x_k))
+
+    where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th
+    root of ``He_n``, and then scaling the results to get the right value
+    when integrating 1.
+
+    """
+    ideg = int(deg)
+    if ideg != deg or ideg < 1:
+        raise ValueError("deg must be a non-negative integer")
+
+    # first approximation of roots. We use the fact that the companion
+    # matrix is symmetric in this case in order to obtain better zeros.
+    c = np.array([0]*deg + [1])
+    m = hermecompanion(c)
+    x = la.eigvals(m)
+    x.sort()
+
+    # improve roots by one application of Newton
+    dy = hermeval(x, c)
+    df = hermeval(x, hermeder(c))
+    x -= dy/df
+
+    # compute the weights. We scale the factor to avoid possible numerical
+    # overflow.
+    fm = hermeval(x, c[1:])
+    fm /= np.abs(fm).max()
+    df /= np.abs(df).max()
+    w = 1/(fm * df)
+
+    # for Hermite_e we can also symmetrize
+    w = (w + w[::-1])/2
+    x = (x - x[::-1])/2
+
+    # scale w to get the right value
+    w *= np.sqrt(2*np.pi) / w.sum()
+
+    return x, w
+
+
+def hermeweight(x):
+    """Weight function of the Hermite_e polynomials.
+
+    The weight function for which the Hermite_e polynomials are orthogonal.
+    In this case the weight function is ``exp(-.5*x**2)``. Note that the
+    Hermite_e polynomials are not normalized.
+
+    Parameters
+    ----------
+    x : array_like
+       Values at which the weight function will be computed.
+
+    Returns
+    -------
+    w : ndarray
+       The weight function at `x`.
+
+    """
+    w = np.exp(-.5*x**2)
+    return w
+
+
 #
 # HermiteE series class
 #
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index 5ab2d23ae..872829fbc 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -39,6 +39,9 @@ Misc Functions
 - `lagvander` -- Vandermonde-like matrix for Laguerre polynomials.
 - `lagvander2d` -- Vandermonde-like matrix for 2D power series.
 - `lagvander3d` -- Vandermonde-like matrix for 3D power series.
+- `laggauss` -- Gauss-Laguerre quadrature, points and weights.
+- `lagweight` -- Laguerre weight function.
+- `lagcompanion` -- symmetrized companion matrix in Laguerre form.
 - `lagfit` -- least-squares fit returning a Laguerre series.
 - `lagtrim` -- trim leading coefficients from a Laguerre series.
 - `lagline` -- Laguerre series of given straight line.
@@ -66,7 +69,8 @@ __all__ = ['lagzero', 'lagone', 'lagx', 'lagdomain', 'lagline',
     'lagadd', 'lagsub', 'lagmulx', 'lagmul', 'lagdiv', 'lagpow',
     'lagval', 'lagder', 'lagint', 'lag2poly', 'poly2lag', 'lagfromroots',
     'lagvander', 'lagfit', 'lagtrim', 'lagroots', 'Laguerre', 'lagval2d',
-    'lagval3d', 'laggrid2d', 'laggrid3d', 'lagvander2d', 'lagvander3d']
+    'lagval3d', 'laggrid2d', 'laggrid3d', 'lagvander2d', 'lagvander3d',
+    'lagcompanion', 'laggauss', 'lagweight']
 
 lagtrim = pu.trimcoef
 
@@ -1504,6 +1508,89 @@ def lagroots(cs):
     return r
 
 
+def laggauss(deg):
+    """Gauss Laguerre quadrature.
+
+    Computes the sample points and weights for Gauss-Laguerre quadrature.
+    These sample points and weights will correctly integrate polynomials of
+    degree ``2*deg - 1`` or less over the interval ``[0, inf]`` with the
+    weight function ``f(x) = exp(-x)``.
+
+    Parameters
+    ----------
+    deg : int
+        Number of sample points and weights. It must be >= 1.
+
+    Returns
+    -------
+    x : ndarray
+        1-D ndarray containing the sample points.
+    y : ndarray
+        1-D ndarray containing the weights.
+
+    Notes
+    -----
+    The results have only been tested up to degree 100. Higher degrees may
+    be problematic. The weights are determined by using the fact that
+
+          w = c / (L'_n(x_k) * L_{n-1}(x_k))
+
+    where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th
+    root of ``L_n``, and then scaling the results to get the right value
+    when integrating 1.
+
+    """
+    ideg = int(deg)
+    if ideg != deg or ideg < 1:
+        raise ValueError("deg must be a non-negative integer")
+
+    # first approximation of roots. We use the fact that the companion
+    # matrix is symmetric in this case in order to obtain better zeros.
+    c = np.array([0]*deg + [1])
+    m = lagcompanion(c)
+    x = la.eigvals(m)
+    x.sort()
+
+    # improve roots by one application of Newton
+    dy = lagval(x, c)
+    df = lagval(x, lagder(c))
+    x -= dy/df
+
+    # compute the weights. We scale the factor to avoid possible numerical
+    # overflow.
+    fm = lagval(x, c[1:])
+    fm /= np.abs(fm).max()
+    df /= np.abs(df).max()
+    w = 1/(fm * df)
+
+    # scale w to get the right value, 1 in this case
+    w /= w.sum()
+
+    return x, w
+
+
+def lagweight(x):
+    """Weight function of the Laguerre polynomials.
+
+    The weight function for which the Laguerre polynomials are orthogonal.
+    In this case the weight function is ``exp(-x)``. Note that the Laguerre
+    polynomials are not normalized, indeed, may be much greater than the
+    normalized versions.
+
+    Parameters
+    ----------
+    x : array_like
+       Values at which the weight function will be computed.
+
+    Returns
+    -------
+    w : ndarray
+       The weight function at `x`.
+
+    """
+    w = np.exp(-x)
+    return w
+
 #
 # Laguerre series class
 #
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index b9a74501c..93ba1d348 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -40,6 +40,9 @@ Misc Functions
 - `legvander` -- Vandermonde-like matrix for Legendre polynomials.
 - `legvander2d` -- Vandermonde-like matrix for 2D power series.
 - `legvander3d` -- Vandermonde-like matrix for 3D power series.
+- `leggauss` -- Gauss-Legendre quadrature, points and weights.
+- `legweight` -- Legendre weight function.
+- `legcompanion` -- symmetrized companion matrix in Legendre form.
 - `legfit` -- least-squares fit returning a Legendre series.
 - `legtrim` -- trim leading coefficients from a Legendre series.
 - `legline` -- Legendre series representing given straight line.
@@ -67,7 +70,8 @@ __all__ = ['legzero', 'legone', 'legx', 'legdomain', 'legline',
     'legadd', 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval',
     'legder', 'legint', 'leg2poly', 'poly2leg', 'legfromroots',
     'legvander', 'legfit', 'legtrim', 'legroots', 'Legendre','legval2d',
-    'legval3d', 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d']
+    'legval3d', 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d',
+    'legcompanion', 'leggauss', 'legweight']
 
 legtrim = pu.trimcoef
 
@@ -1510,6 +1514,92 @@ def legroots(cs):
     return r
 
 
+def leggauss(deg):
+    """Gauss Legendre quadrature.
+
+    Computes the sample points and weights for Gauss-Legendre quadrature.
+    These sample points and weights will correctly integrate polynomials of
+    degree ``2*deg - 1`` or less over the interval ``[-1, 1]`` with the
+    weight function ``f(x) = 1``.
+
+    Parameters
+    ----------
+    deg : int
+        Number of sample points and weights. It must be >= 1.
+
+    Returns
+    -------
+    x : ndarray
+        1-D ndarray containing the sample points.
+    y : ndarray
+        1-D ndarray containing the weights.
+
+    Notes
+    -----
+    The results have only been tested up to degree 100. Higher degrees may
+    be problematic. The weights are determined by using the fact that
+
+          w = c / (L'_n(x_k) * L_{n-1}(x_k))
+
+    where ``c`` is a constant independent of ``k`` and ``x_k`` is the k'th
+    root of ``L_n``, and then scaling the results to get the right value
+    when integrating 1.
+
+    """
+    ideg = int(deg)
+    if ideg != deg or ideg < 1:
+        raise ValueError("deg must be a non-negative integer")
+
+    # first approximation of roots. We use the fact that the companion
+    # matrix is symmetric in this case in order to obtain better zeros.
+    c = np.array([0]*deg + [1])
+    m = legcompanion(c)
+    x = la.eigvals(m)
+    x.sort()
+
+    # improve roots by one application of Newton
+    dy = legval(x, c)
+    df = legval(x, legder(c))
+    x -= dy/df
+
+    # compute the weights. We scale the factor to avoid possible numerical
+    # overflow.
+    fm = legval(x, c[1:])
+    fm /= np.abs(fm).max()
+    df /= np.abs(df).max()
+    w = 1/(fm * df)
+
+    # for Legendre we can also symmetrize
+    w = (w + w[::-1])/2
+    x = (x - x[::-1])/2
+
+    # scale w to get the right value
+    w *= 2. / w.sum()
+
+    return x, w
+
+
+def legweight(x):
+    """Weight function of the Legendre polynomials.
+
+    The weight function for which the Legendre polynomials are orthogonal.
+    In this case the weight function is simply one. Note that the Legendre
+    polynomials are not normalized.
+
+    Parameters
+    ----------
+    x : array_like
+       Values at which the weight function will be computed.
+
+    Returns
+    -------
+    w : ndarray
+       The weight function at `x`.
+
+    """
+    w = x*0.0 + 1.0
+    return w
+
 #
 # Legendre series class
 #