program montecarlo
  implicit none
  integer,parameter :: kosuu=1024 ! 発生させる乱数の数
  real :: p,x,y,s,ysouwa,y2souwa,sekibun,bunsan,sdom
  real,parameter :: x_kaishi=0.0,x_shuryo=1.0 ! 積分区間の始点と終点
  real,parameter :: x_haba=x_shuryo-x_kaishi
  integer :: i
  ysouwa = 0.0
  y2souwa = 0.0
  do i = 1,kosuu
     call random_number(p)  ! 0から1の範囲でランダムに実数 p を発生させ
     x = x_kaishi+p*x_haba  ! p に基づいて x_kaishi から x_shuryo までの実数 x を作る
     y = 1.0/(1.0+x**2)     ! x における被積分関数の値を計算
     s = y*x_haba           ! 高さ y、幅 x_haba の長方形として定積分を近似する
     ysouwa = ysouwa+s      !「長方形の面積」の総和をとる
     y2souwa = y2souwa+s**2 !「長方形の面積」の2乗の総和をとる
  end do
  sekibun = ysouwa/real(kosuu) ! 「長方形の面積」の「期待値」として定積分を求める
  bunsan = y2souwa/real(kosuu)-sekibun**2 ! 「長方形の面積」の不偏でない「分散」
  sdom = sqrt(bunsan/real(kosuu-1))       ! 「長方形の面積」の「平均値の標準偏差」
  write(*,*)kosuu,sekibun,sekibun*4.0,abs(sekibun-atan(1.0)),sdom
end program montecarlo