program main

!*****************************************************************************80
!
!! MAIN is the main program for MC.
!
!  Discussion:
!
!    Things you can vary include:
!
!    * the function to be integrated (see the list below)
!    * the number of samples used, SAMPLE_NUM.
!    * the spatial dimension D_NUM (try 1, 2, 4)
!
!    The function to be integrated is:
!
!      F = Product ( 1 <= DIM <= DIM_NUM ) ( 2 * abs ( 2 * X(DIM) - 1 ) ) 
!      Exact integral is 1.
!
!  You will need the "USE" statement here in order to access MPI functions
!  and variables.
   
  integer ( kind = 4 ), parameter :: dim_num = 4

  real    ( kind = 8 ) f
  real    ( kind = 8 ) q
  real    ( kind = 8 ) q_error
  real    ( kind = 8 ), parameter :: q_exact = 1.0D+00
  integer ( kind = 4 ) sample
  integer ( kind = 4 ), parameter :: sample_num = 1000
  integer ( kind = 4 ) seed
  real    ( kind = 8 ) x(dim_num)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'MC'
  write ( *, '(a)' ) '  FORTRAN90 version.'

  seed = 123456789

  q = 0.0D+00

  do sample = 1, sample_num

    call r8vec_uniform_01 ( dim_num, seed, x )

    q = q + f ( dim_num, x )

  end do

  q = q / real ( sample_num, kind = 8 )
  q_error = abs ( q - q_exact )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   Samples Dimension      Estimate          Error'
  write ( *, '(a)' ) ' '

  write ( *, '(2x,i8,2x,i8,2x,f16.10,2x,g16.6)' ) &
    sample_num, dim_num, q, q_error

  stop
end
function f ( dim_num, x )

!*****************************************************************************80
!
!! F evaluates the function F(X) which we are integrating.
!
!  Modified:
!
!    04 September 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the number of spatial dimensions.
!
!    Input, real ( kind = 8 ) X(DIM_NUM), the point at which to evaluate F.
!
!    Output, real ( kind = 8 ) F, the value of F(X).
!
  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) f
  real ( kind = 8 ) x(dim_num)

  f = product ( 2.0D+00 * abs ( 2.0D+00 * x(1:dim_num) - 1.0D+00 ) )

  return
end
subroutine r8vec_uniform_01 ( n, seed, r )

!*****************************************************************************80
!
!! R8VEC_UNIFORM_01 returns a unit pseudorandom vector.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the vector.
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which 
!    should NOT be 0.  On output, SEED has been updated.
!
!    Output, real ( kind = 8 ) R(N), the vector of pseudorandom values.
!
  implicit none

  integer ( kind = 4 ) n

  integer ( kind = 4 ) i
  integer ( kind = 4 ), parameter :: i4_huge = 2147483647
  integer ( kind = 4 ) k
  integer ( kind = 4 ) seed
  real    ( kind = 8 ) r(n)

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop
  end if

  do i = 1, n

    k = seed / 127773

    seed = 16807 * ( seed - k * 127773 ) - k * 2836

    if ( seed < 0 ) then
      seed = seed + i4_huge
    end if

    r(i) = real ( seed, kind = 8 ) * 4.656612875D-10

  end do

  return
end