legzo Subroutine

subroutine legzo ( n, x, w )

  !*****************************************************************************80
  !
  !! LEGZO computes the zeros of Legendre polynomials, and integration weights.
  !
  !  Discussion:
  !
  !    This procedure computes the zeros of Legendre polynomial Pn(x) in the 
  !    interval [-1,1], and the corresponding weighting coefficients for 
  !    Gauss-Legendre 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:
  !
  !    13 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 polynomial.
  !
  !    Output, real ( kind = 8 ) X(N), W(N), the zeros of the polynomial,
  !    and the corresponding weights.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) fd
  real ( kind = 8 ) gd
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) n0
  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

  n0 = ( n + 1 ) / 2

  do nr = 1, n0

     z = cos ( 3.1415926D+00 * ( nr - 0.25D+00 ) / n )

     do

        z0 = z
        p = 1.0D+00
        do i = 1, nr - 1
           p = p * ( z - x(i))
        end do
        f0 = 1.0D+00
        if ( nr == n0 .and. n /= 2 * int ( n / 2 ) ) then
           z = 0.0D+00
        end if
        f1 = z
        do k = 2, n
           pf = ( 2.0D+00 - 1.0D+00 / k ) * z * f1 &
                - ( 1.0D+00 - 1.0D+00 / k ) * f0
           pd = k * ( f1 - z * pf ) / ( 1.0D+00 - z * z )
           f0 = f1
           f1 = pf
        end do

        if ( z == 0.0D+00 ) then
           exit
        end if

        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 ( abs ( z - z0 ) < abs ( z ) * 1.0D-15 ) then
           exit
        end if

     end do

     x(nr) = z
     x(n+1-nr) = - z
     w(nr) = 2.0D+00 / ( ( 1.0D+00 - z * z ) * pd * pd )
     w(n+1-nr) = w(nr)

  end do

  return
end subroutine legzo