Monte Carlo で円周率

半径１の１／４円中に入るか入らないかで円周率を求める古典的な例です。

一千万回の試行をしているけれども、組み込み乱数がちゃんとしているかは不明ｗ
試行回数Nに対して、おおむね 1/sqrt(N) に比例して誤差が減るはずなので目安の数字も出力しています。

ソース・プログラム

program MonteCarlo
  implicit none
  real, parameter :: pi = 4.0 * atan(1.0)
  integer, parameter :: nn = 10**7
  real :: x(nn), y(nn), pp
  call random_seed() 
  call random_number(x)
  call random_number(y)
  pp = 4 * count( x**2 + y**2 < 1.0 ) / real(nn)
  print '(3f12.7)', pp, abs(pp - pi),  1.0 / sqrt( real(nn) )
  stop 
end program MonteCarlo