function [ xtab, weight ] = gen_laguerre_compute ( order, alpha )

%% GEN_LAGUERRE_COMPUTE computes a generalized Gauss-Laguerre quadrature rule.
%
%  Discussion:
%
%    The integration interval is [ 0, +oo ).
%
%    The weight function w(x) = EXP ( - X ) * X**ALPHA.
%
%
%    If the integral to approximate is:
%
%        Integral ( 0 <= X < +oo ) EXP ( - X ) * X**ALPHA * F(X) dX
%
%    then the quadrature rule is:
%
%      Sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( XTAB(I) )
%
%
%    If the integral to approximate is:
%
%        Integral ( 0 <= X < +oo ) X**ALPHA * F(X) dX
%
%    then the quadrature rule is:
%
%      Sum ( 1 <= I <= ORDER ) WEIGHT(I) * EXP(XTAB(I)) * F ( XTAB(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2008
%
%  Author:
%
%    FORTRAN77 original version by Arthur Stroud, Don Secrest
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the quadrature rule to be computed.
%    ORDER must be at least 1.
%
%    Input, real ALPHA, the exponent of the X factor.
%    Set ALPHA = 0.0 for the simplest rule.
%    ALPHA must be nonnegative.
%
%    Output, real XTAB(ORDER), the Gauss-Laguerre abscissas.
%
%    Output, real WEIGHT(ORDER), the Gauss-Laguerre weights.
%

%
%  Set the recursion coefficients.
%
  for i = 1 : order
    b(i) = ( alpha + 2 * i - 1 );
  end

  for i = 1 : order
    c(i) = ( i - 1 ) * ( alpha + i - 1 );
  end

  cc = r8_gamma ( alpha + 1.0 ) * prod ( c(2:order) );

  for i = 1 : order
%
%  Compute an estimate for the root.
%
    if ( i == 1 )

      x = ( 1.0 + alpha ) * ( 3.0 + 0.92 * alpha ) ...
        / ( 1.0 + 2.4 * order + 1.8 * alpha );

    elseif ( i == 2 )

      x = x + ( 15.0 + 6.25 * alpha ) ...
        / ( 1.0 + 0.9 * alpha + 2.5 * order );

    else

      r1 = ( 1.0 + 2.55 * ( i - 2 ) ) / ( 1.9 * ( i - 2 ) );

      r2 = 1.26 * ( i - 2 ) * alpha / ( 1.0 + 3.5 * ( i - 2 ) );

      ratio = ( r1 + r2 ) / ( 1.0 + 0.3 * alpha );

      x = x + ratio * ( x - xtab(i-2) );

    end
%
%  Use iteration to find the root.
%
    [ x, dp2, p1 ] = gen_laguerre_root ( x, order, alpha, b, c );
%
%  Set the abscissa and weight.
%
    xtab(i) = x;
    weight(i) = ( cc / dp2 ) / p1;

  end