program sphere_volume implicit none integer, parameter :: dp = kind(0.0d0) real(dp), parameter :: pi = acos(-1.0_dp) integer :: n, n_dim, n_runs, n_samples real(dp) :: volume, mean, std_dev real(dp), allocatable :: inside_frac(:) character(len=20) :: arg if (command_argument_count() /= 3) then stop "Usage: ./sphere_volume n_dim n_runs n_samples" end if call get_command_argument(1, arg) read(arg, *) n_dim call get_command_argument(2, arg) read(arg, *) n_runs call get_command_argument(3, arg) read(arg, *) n_samples call set_random_seed() allocate(inside_frac(n_runs)) do n = 1, n_runs call run_monte_carlo(n_dim, n_samples, inside_frac(n)) end do mean = sum(inside_frac) / n_runs std_dev = sqrt(sum((inside_frac - mean)**2) / n_runs) mean = mean * 2.0_dp**n_dim std_dev = std_dev * 2.0_dp**n_dim print *, n_samples, mean, std_dev, solution(1.0_dp, n_dim) contains subroutine set_random_seed() integer, allocatable :: seed(:) integer :: time, i, size_seed call random_seed(size=size_seed) allocate(seed(size_seed)) call system_clock(time) do i = 1, size_seed seed(i) = ishftc(time, i*8) end do call random_seed(put=seed) end subroutine set_random_seed subroutine run_monte_carlo(n_dim, n_samples, inside_frac) integer, intent(in) :: n_dim, n_samples real(dp), intent(out) :: inside_frac real(dp) :: x(n_dim) integer :: n, inside_count inside_count = 0 do n = 1, n_samples call random_number(x) if (sum(x**2) < 1.0_dp) then inside_count = inside_count + 1 end if end do inside_frac = inside_count / real(n_samples, dp) end subroutine run_monte_carlo function solution(radius, ndim) result(volume) real(dp), intent(in) :: radius integer, intent(in) :: ndim real(dp) :: volume volume = radius**ndim * pi**(0.5_dp * ndim) / & gamma(0.5_dp * ndim + 1.0_dp) end function solution end program