1program answer
2 implicit none
3
4 integer, parameter :: max_memory_mb = 512 ! 最大メモリ使用量 [MB]
5 real(8), parameter :: pi = atan(1.0_8) * 4
6
7 integer :: ndims, trial
8 real(8) :: approx, exact
9
10 ! 次元数と試行回数を求める
11 read(*, *) ndims, trial
12
13 call random_seed_clock()
14 approx = volume(ndims, trial) * 2**ndims
15 exact = pi**(ndims * 0.5_8) / gamma(ndims * 0.5_8 + 1)
16
17 write(*, fmt='("approximation = ", e20.8)') approx
18 write(*, fmt='("exact value = ", e20.8)') exact
19 write(*, fmt='("relative error = ", e20.8)') abs(1 - approx / exact)
20
21 stop
22contains
23
24 ! 超球の体積を求める
25 real(8) function volume(ndims, trial)
26 implicit none
27 integer, intent(in) :: ndims
28 integer, intent(in) :: trial
29
30 integer :: i, j, k, niter, count, nmax
31 real(8), allocatable :: x(:, :)
32
33 ! 乱数生成個数の最大値 (メモリ使用量が大きくなり過ぎないように)
34 nmax = max_memory_mb * 2**20 / (8 * ndims)
35
36 ! 反復回数
37 niter = trial / nmax
38
39 ! モンテカルロ積分
40 allocate(x(ndims, nmax))
41
42 count = 0
43 do k = 1, niter
44 call random_number(x)
45 count = count + count_inside(x)
46 enddo
47
48 call random_number(x)
49 count = count + count_inside(x(:, 1:mod(trial, nmax)))
50
51 volume = real(count, kind=8) / real(trial, kind=8)
52
53 deallocate(x)
54
55 endfunction volume
56
57 ! 試行
58 integer function count_inside(x)
59 real(8), intent(in) :: x(:, :)
60
61 integer :: i, j, n, m
62 real(8) :: r
63
64 n = size(x, 1)
65 m = size(x, 2)
66
67 ! 条件を満たす要素を数える
68 count_inside = 0
69 do j = 1, m
70 r = sum(x(:, j)**2)
71
72 if(r < 1.0_8) then
73 count_inside = count_inside + 1
74 endif
75 enddo
76
77 endfunction count_inside
78
79 ! 乱数のseedをシステムクロックに応じて変更
80 subroutine random_seed_clock()
81 implicit none
82 integer :: nseed, clock
83 integer, allocatable :: seed(:)
84
85 call system_clock(clock)
86
87 call random_seed(size=nseed)
88 allocate(seed(nseed))
89
90 seed = clock
91 call random_seed(put=seed)
92
93 deallocate(seed)
94 endsubroutine random_seed_clock
95
96endprogram answer