function table = laguerre_linear_product ( p, e )

%*****************************************************************************80
%
%% LAGUERRE_LINEAR_PRODUCT computes a linearly weighted Laguerre product table.
%
%  Discussion:
%
%    Let L(i)(X) represent the Laguerre polynomial of degree i.  
%
%    For polynomial chaos applications, it is of interest to know the
%    value of the integrals of products of X with every possible pair
%    of basis functions.  That is, we'd like to form
%
%      Tij = Integral ( 0 <= X < +oo ) X^E * L(i)(X) * L(j)(X) * exp (-X) dx
%
%  Modified:
%
%    09 May 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the maximum degree of the polyonomial factors.
%    0 <= P.
%
%    Input, integer E, the exponent of X in the integrand.
%    0 <= E.
%
%    Output, real TABLE(P+1,P+1), the table of integrals.  TABLE(I,J)
%    represents the weighted integral of X^E * L(i+1)(X) * L(j+1)(X).
%
  table(1:p+1,1:p+1) = 0.0;

  order = p + 1 + floor ( ( e + 1 ) / 2 );

  [ x_table, w_table ] = laguerre_quadrature_rule ( order );

  for k = 1 : order

    x = x_table(k);
    l_table = laguerre_polynomial ( p, x );
%
%  The following formula is an outer product in L_TABLE.
%
    if ( e == 0 )
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) *       l_table(1:p+1)' * l_table(1:p+1);
    else
      table(1:p+1,1:p+1) = table(1:p+1,1:p+1) ...
        + w_table(k) * x^e * l_table(1:p+1)' * l_table(1:p+1);
    end

  end

  return
end

