lagzo Subroutine

subroutine lagzo ( n, x, w )

  !*****************************************************************************80
  !
  !! LAGZO computes zeros of the Laguerre polynomial, and integration weights.
  !
  !  Discussion:
  !
  !    This procedure computes the zeros of Laguerre polynomial Ln(x) in the 
  !    interval [0,∞], and the corresponding weighting coefficients for 
  !    Gauss-Laguerre integration.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    07 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) N, the order of the Laguerre polynomial.
  !
  !    Output, real ( kind = 8 ) X(N), the zeros of the Laguerre polynomial.
  !
  !    Output, real ( kind = 8 ) W(N), the weighting coefficients.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) fd
  real ( kind = 8 ) gd
  real ( kind = 8 ) hn
  integer ( kind = 4 ) i
  integer ( kind = 4 ) it
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) nr
  real ( kind = 8 ) p
  real ( kind = 8 ) pd
  real ( kind = 8 ) pf
  real ( kind = 8 ) q
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) wp
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) z
  real ( kind = 8 ) z0

  hn = 1.0D+00 / real ( n, kind = 8 )

  do nr = 1, n

     if ( nr == 1 ) then
        z = hn
     else
        z = x(nr-1) + hn * nr ** 1.27D+00
     end if

     it = 0

     do

        it = it + 1
        z0 = z
        p = 1.0D+00
        do i = 1, nr - 1
           p = p * ( z - x(i) )
        end do

        f0 = 1.0D+00
        f1 = 1.0D+00 - z
        do k = 2, n
           pf = (( 2.0D+00 * k - 1.0D+00 - z ) * f1 &
                - ( k - 1.0D+00 ) * f0 ) / k
           pd = k / z * ( pf - f1 )
           f0 = f1
           f1 = pf
        end do

        fd = pf / p

        q = 0.0D+00
        do i = 1, nr - 1
           wp = 1.0D+00
           do j = 1, nr - 1
              if ( j /= i ) then
                 wp = wp * ( z - x(j) )
              end if
           end do
           q = q + wp
        end do

        gd = ( pd - q * fd ) / p
        z = z - fd / gd

        if ( 40 < it .or. abs ( ( z - z0 ) / z ) <= 1.0D-15 ) then
           exit
        end if

     end do

     x(nr) = z
     w(nr) = 1.0D+00 / ( z * pd * pd )

  end do

  return
end subroutine lagzo