subroutine chebyshev1_exactness ( n, x, w, p_max )

!*****************************************************************************80
!
!! chebyshev1_exactness(): monomial exactness for the Chebyshev1 integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    27 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of points in the rule.
!
!    real ( kind = rk ) X(N), the quadrature points.
!
!    real ( kind = rk ) W(N), the quadrature weights.
!
!    integer P_MAX, the maximum exponent.
!    0 <= P_MAX.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) e
  integer p
  integer p_max
  real ( kind = rk ) q
  real ( kind = rk ) s
  real ( kind = rk ) v(n)
  real ( kind = rk ) w(n)
  real ( kind = rk ) x(n)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Quadrature rule for the Chebyshev1 integral.'
  write ( *, '(a,i3)' ) '  Rule has order N = ', n
  write ( *, '(a)' ) '  Degree          Relative Error'
  write ( *, '(a)' ) ''

  do p = 0, p_max

    call chebyshev1_integral ( p, s )

    v(1:n) = x(1:n) ** p

    q = dot_product ( w, v )

    if ( s == 0.0D+00 ) then
      e = abs ( q )
    else
      e = abs ( q - s ) / abs ( s )
    end if

    write ( *, '(2x,i6,2x,f24.16)' ) p, e

  end do

  return
end
subroutine chebyshev1_integral ( expon, exact )

!*****************************************************************************80
!
!! CHEBYSHEV1_INTEGRAL evaluates the Chebyshev type 1 integral of a monomial.
!
!  Discussion:
!
!    The integral:
!
!      integral ( -1 <= x <= +1 ) x^n / sqrt ( 1 - x * x ) dx
!
!    This routine is given the value of the exponent, and returns the
!    exact value of the integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    22 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer EXPON, the exponent.
!
!  Output:
!
!    real ( kind = rk ) EXACT, the value of the exact integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) bot
  real ( kind = rk ) exact
  integer expon
  integer i
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) top
!
!  Get the exact value of the integral.
!
  if ( mod ( expon, 2 ) == 0 ) then

    top = 1
    bot = 1
    do i = 2, expon, 2
      top = top * ( i - 1 )
      bot = bot *   i
    end do
    
    exact = r8_pi * real ( top, kind = rk ) / real ( bot, kind = rk )

  else

    exact = 0.0D+00
    
  end if

  return
end
subroutine chebyshev2_exactness ( n, x, w, p_max )

!*****************************************************************************80
!
!! CHEBYSHEV2_EXACTNESS: monomial exactness for the Chebyshev2 integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    28 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of points in the rule.
!
!    real ( kind = rk ) X(N), the quadrature points.
!
!    real ( kind = rk ) W(N), the quadrature weights.
!
!    integer P_MAX, the maximum exponent.
!    0 <= P_MAX.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) e
  integer p
  integer p_max
  real ( kind = rk ) q
  real ( kind = rk ) s
  real ( kind = rk ) v(n)
  real ( kind = rk ) w(n)
  real ( kind = rk ) x(n)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Quadrature rule for the Chebyshev2 integral.'
  write ( *, '(a,i3)' ) '  Rule has order N = ', n
  write ( *, '(a)' ) '  Degree          Relative Error'
  write ( *, '(a)' ) ''

  do p = 0, p_max

    call chebyshev2_integral ( p, s )

    v(1:n) = x(1:n) ** p

    q = dot_product ( w, v )

    if ( s == 0.0D+00 ) then
      e = abs ( q )
    else
      e = abs ( q - s ) / abs ( s )
    end if

    write ( *, '(2x,i6,2x,f24.16)' ) p, e

  end do

  return
end
subroutine chebyshev2_integral ( expon, exact )

!*****************************************************************************80
!
!! CHEBYSHEV2_INTEGRAL evaluates a monomial Chebyshev type 2 integral.
!
!  Discussion:
!
!    The integral:
!
!      integral ( -1 <= x <= +1 ) x^n * sqrt ( 1 - x * x ) dx
!
!    This routine is given the value of the exponent, and returns the
!    exact value of the integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    22 February 2008
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer EXPON, the exponent.
!
!  Output:
!
!    real ( kind = rk ) EXACT, the value of the exact integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) bot
  real ( kind = rk ) exact
  integer expon
  integer i
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) top
!
!  Get the exact value of the integral.
!
  if ( mod ( expon, 2 ) == 0 ) then

    top = 1
    bot = 1
    do i = 2, expon, 2
      top = top * ( i - 1 )
      bot = bot *   i
    end do

    bot = bot * real ( expon + 2, kind = rk )

    exact = r8_pi * real ( top, kind = rk ) / real ( bot, kind = rk )

  else

    exact = 0.0D+00
    
  end if

  return
end
subroutine gegenbauer_exactness ( n, x, w, p_max, lambda )

!*****************************************************************************80
!
!! GEGENBAUER_EXACTNESS: monomial exactness for the Gegenbauer integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    13 January 2016
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of points in the rule.
!
!    real ( kind = rk ) X(N), the quadrature points.
!
!    real ( kind = rk ) W(N), the quadrature weights.
!
!    integer P_MAX, the maximum exponent.
!    0 <= P_MAX.
!
!    real ( kind = rk ) LAMBDA, the parameter.
!    -1/2 < LAMBDA,
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) e
  real ( kind = rk ) lambda
  integer p
  integer p_max
  real ( kind = rk ) q
  real ( kind = rk ) s
  real ( kind = rk ) v(n)
  real ( kind = rk ) w(n)
  real ( kind = rk ) x(n)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Quadrature rule for the Gegenbauer integral.'
  write ( *, '(a,g14.6)' ) '  Lambda = ', lambda
  write ( *, '(a,i3)' ) '  Rule has order N = ', n
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Degree          Relative Error'
  write ( *, '(a)' ) ''

  do p = 0, p_max

    call gegenbauer_integral ( p, lambda, s )

    v(1:n) = x(1:n) ** p

    q = dot_product ( w, v )

    if ( s == 0.0D+00 ) then
      e = abs ( q )
    else
      e = abs ( q - s ) / abs ( s )
    end if

    write ( *, '(2x,i6,2x,f24.16)' ) p, e

  end do

  return
end
subroutine gegenbauer_integral ( p, lambda, s )

!*****************************************************************************80
!
!! GEGENBAUER_INTEGRAL evaluates a monomial integral with Gegenbauer weight.
!
!  Discussion:
!
!    The integral:
!
!      integral ( -1 <= x < +1 ) x^p * ( 1 - x^2 )^(lambda-1/2) dx
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    10 January 2016
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer P, the exponent.
!    0 <= P.
!
!    real ( kind = rk ) LAMBDA, the exponent term.
!    -1/2 < LAMBDA.
!
!  Output:
!
!    real ( kind = rk ) S, the value of the integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) lambda
  integer p
  real ( kind = rk ) s

  if ( mod ( p, 2 ) == 0 ) then
    s = gamma ( p / 2.0D+00 + 0.5D+00 ) * gamma ( lambda + 0.5D+00 ) &
      / gamma ( p / 2.0D+00 + lambda + 1.0D+00 )
  else
    s = 0.0D+00
  end if

  return
end
subroutine hermite_exactness ( n, x, w, p_max )

!*****************************************************************************80
!
!! HERMITE_EXACTNESS: monomial exactness for the Hermite integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    16 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of points in the rule.
!
!    real ( kind = rk ) X(N), the quadrature points.
!
!    real ( kind = rk ) W(N), the quadrature weights.
!
!    integer P_MAX, the maximum exponent.
!    0 <= P_MAX.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) e
  integer p
  integer p_max
  real ( kind = rk ) q
  real ( kind = rk ) s
  real ( kind = rk ) v(n)
  real ( kind = rk ) w(n)
  real ( kind = rk ) x(n)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Quadrature rule for the Hermite integral.'
  write ( *, '(a,i3)' ) '  Rule has order N = ', n
  write ( *, '(a)' ) '  Degree          Relative Error'
  write ( *, '(a)' ) ''

  do p = 0, p_max

    call hermite_integral ( p, s )

    v(1:n) = x(1:n) ** p

    q = dot_product ( w, v )

    if ( s == 0.0D+00 ) then
      e = abs ( q )
    else
      e = abs ( q - s ) / abs ( s )
    end if

    write ( *, '(2x,i6,2x,f24.16)' ) p, e

  end do

  return
end
subroutine hermite_integral ( p, s )

!*****************************************************************************80
!
!! HERMITE_INTEGRAL evaluates a monomial Hermite integral.
!
!  Discussion:
!
!    Integral ( -oo < x < +oo ) x^p exp(-x^2) dx
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    18 November 2011
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer P, the exponent. 
!    0 <= P.
!
!  Output:
!
!    real ( kind = rk ) S, the value of the integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer p
  real ( kind = rk ) r8_factorial2
  real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00
  real ( kind = rk ) s

  if ( mod ( p, 2 ) == 0 ) then
    s = r8_factorial2 ( p - 1 ) * sqrt ( r8_pi ) / 2.0D+00 ** ( p / 2 )
  else
    s = 0.0D+00
  end if

  return
end
subroutine laguerre_exactness ( n, x, w, p_max )

!*****************************************************************************80
!
!! LAGUERRE_EXACTNESS: monomial exactness for the Laguerre integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    16 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of points in the rule.
!
!    real ( kind = rk ) X(N), the quadrature points.
!
!    real ( kind = rk ) W(N), the quadrature weights.
!
!    integer P_MAX, the maximum exponent.
!    0 <= P_MAX.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) e
  integer p
  integer p_max
  real ( kind = rk ) q
  real ( kind = rk ) s
  real ( kind = rk ) v(n)
  real ( kind = rk ) w(n)
  real ( kind = rk ) x(n)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Quadrature rule for the Laguerre integral.'
  write ( *, '(a,i3)' ) '  Rule has order N = ', n
  write ( *, '(a)' ) '  Degree          Relative Error'
  write ( *, '(a)' ) ''

  do p = 0, p_max

    call laguerre_integral ( p, s )

    v(1:n) = x(1:n) ** p

    q = dot_product ( w, v )

    if ( s == 0.0D+00 ) then
      e = abs ( q )
    else
      e = abs ( q - s ) / abs ( s )
    end if

    write ( *, '(2x,i6,2x,f24.16)' ) p, e

  end do

  return
end
subroutine laguerre_integral ( p, s )

!*****************************************************************************80
!
!! LAGUERRE_INTEGRAL evaluates a monomial Laguerre integral.
!
!  Discussion:
!
!    The integral:
!
!      integral ( 0 <= x < +oo ) x^p * exp ( -x ) dx
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 March 2012
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer P, the exponent.
!    0 <= P.
!
!  Output:
!
!    real ( kind = rk ) S, the value of the integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer p
  real ( kind = rk ) r8_factorial
  real ( kind = rk ) s

  s = r8_factorial ( p )

  return
end
subroutine legendre_exactness ( n, x, w, p_max )

!*****************************************************************************80
!
!! LEGENDRE_EXACTNESS: monomial exactness for the Legendre integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    16 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of points in the rule.
!
!    real ( kind = rk ) X(N), the quadrature points.
!
!    real ( kind = rk ) W(N), the quadrature weights.
!
!    integer P_MAX, the maximum exponent.
!    0 <= P_MAX.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) e
  integer p
  integer p_max
  real ( kind = rk ) q
  real ( kind = rk ) s
  real ( kind = rk ) v(n)
  real ( kind = rk ) w(n)
  real ( kind = rk ) x(n)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '  Quadrature rule for the Legendre integral.'
  write ( *, '(a,i3)' ) '  Rule has order N = ', n
  write ( *, '(a)' ) '  Degree          Relative Error'
  write ( *, '(a)' ) ''

  do p = 0, p_max

    call legendre_integral ( p, s )

    v(1:n) = x(1:n) ** p

    q = dot_product ( w, v )

    if ( s == 0.0D+00 ) then
      e = abs ( q )
    else
      e = abs ( q - s ) / abs ( s )
    end if

    write ( *, '(2x,i6,2x,f24.16)' ) p, e

  end do

  return
end
subroutine legendre_integral ( p, s )

!*****************************************************************************80
!
!! LEGENDRE_INTEGRAL evaluates a monomial Legendre integral.
!
!  Discussion:
!
!    To test a Legendre quadrature rule, we use it to approximate the
!    integral of a monomial:
!
!      integral ( -1 <= x <= +1 ) x^p dx
!
!    This routine is given the value of the exponent, and returns the
!    exact value of the integral.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    16 May 2014
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer P, the power.
!
!  Output:
!
!    real ( kind = rk ) S, the value of the exact integral.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer p
  real ( kind = rk ) s
!
!  Get the exact value of the integral.
!
  if ( mod ( p, 2 ) == 0 ) then
  
    s = 2.0D+00 / real ( p + 1, kind = rk )
	
  else
  
    s = 0.0D+00
	
  end if

  return
end
function r8_factorial ( n )

!*****************************************************************************80
!
!! R8_FACTORIAL computes the factorial of N.
!
!  Discussion:
!
!    factorial ( N ) = product ( 1 <= I <= N ) I
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the argument of the factorial function.
!    If N is less than 1, the function value is returned as 1.
!
!  Output:
!
!    real ( kind = rk ) R8_FACTORIAL, the factorial of N.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) r8_factorial
  integer i
  integer n

  r8_factorial = 1.0D+00

  do i = 1, n
    r8_factorial = r8_factorial * real ( i, kind = rk )
  end do

  return
end
function r8_factorial2 ( n )

!*****************************************************************************80
!
!! R8_FACTORIAL2 computes the double factorial function.
!
!  Discussion:
!
!    FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 )  (N even)
!                    = Product ( N * (N-2) * (N-4) * ... * 1 )  (N odd)
!
!  Example:
!
!     N    N!!
!
!     0     1
!     1     1
!     2     2
!     3     3
!     4     8
!     5    15
!     6    48
!     7   105
!     8   384
!     9   945
!    10  3840
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    02 September 2007
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the argument of the double factorial 
!    function.  If N is less than 1, R8_FACTORIAL2 is returned as 1.0.
!
!  Output:
!
!    real ( kind = rk ) R8_FACTORIAL2, the value of N!!.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n
  real ( kind = rk ) r8_factorial2
  real ( kind = rk ) r8_n
  real ( kind = rk ) value

  if ( n < 1 ) then

    value = 1.0D+00

  else

    r8_n = real ( n, kind = rk )
    value = 1.0D+00

    do while ( 1.0D+00 < r8_n )
      value = value * r8_n
      r8_n = r8_n - 2.0D+00
    end do

  end if

  r8_factorial2 = value

  return
end