モンテカルロ法で楕円の面積を求めてみる

なんとなーく CPU をぶん回してみようかなとしょうもないネタを。

ほぼ、マニュアル https://t.co/SpnKeEDZqw を眺めながら取り敢えず動くやつ。
4コアやから、マルチスレッドにせんと意味が無いねぇ。表示ルーチンも毎回モジュロ取るよりタイマー割り込みの方が軽いやろうし。F90 ぐらいの固定形式？ pic.twitter.com/6UiJLt6U40

— 1円切手 (&o53) (@jp_yen) 2018年4月10日

この辺を少し書き直してみました。

MonteCarlo

! 楕円の面積をモンテカルロ法で求める

! 原点を中心とする横長の楕円を考える、第一象限の面積を求め 4倍する
! 長半径, 短半径 Ra, Rb
! 焦点の座標 (F,0), (-F,0)
!
! Ra, Rb ; 長半径, 短半径
! F      ; 中心から楕円の焦点までの距離 Fa(-F, 0), Fb(F, 0)
! Px, Py ; ランダムなポイント P(Px, Py)
! I      ; カウンタ
! Imax   ; 最大試行回数
! J      ; 楕円の範囲内であった数
! K      ; 試行回数
! L      ; 焦点 Fa, Fb から円周上点P までの距離の合計 2 * Ra
! Lab    ; 焦点 Fa, Fb からランダムな点P までの距離の合計
!
! (0,0) - (Ra,Rb) の正方形の中にランダムに点を打ち、楕円の中に
! 入っているか判定する
!
! 初期化
      PROGRAM EllipseArea_MonteCarlo
      IMPLICIT NONE
      INTEGER I, Istep, Imax, J, K
      INTEGER SeedSize
      INTEGER,ALLOCATABLE,DIMENSION(:) :: Seed
      REAL*16 F, Ra, Rb, L, Lab, Px, Py, RaRb4, AT_1, CORRECT

      INTEGER Time_Start, Time_End, t_rate, t_max
      REAL Time_Diff
      REAL CpuTime_Start, CpuTime_End

      ! 実行時間を測る
      CALL SYSTEM_CLOCK(Time_Start, t_rate, t_max)! 開始時刻
      CALL CPU_TIME(CpuTime_Start)
      ! 変数初期化
      Ra = 160
      Rb =  90
      L = 2 * Ra
      J = 0
      F = SQRT( Ra ** 2 - Rb ** 2 )     ! 焦点を求める
      RaRb4   = Ra * Rb * 4.0           ! 4 * Ra * Rb 計算省略用
      AT_1    = ATAN(1.0)               ! arctan(1)   計算省略用
      CORRECT = RaRb4 * AT_1            ! 面積の真値  計算省略用

      Imax = HUGE(Imax)                 ! Int の最大値まで頑張る
      Imax = 50000000                   ! いや、頑張れませんでした。
      Istep= 1000000                    ! 途中経過を表示する単位
      Imax = Istep * INT(Imax/Istep)    ! Imax を突き抜けないように計算しなおし

      ! 乱数初期化
      !
      CALL RANDOM_SEED(SIZE = SeedSize) ! seed のサイズを求める
      WRITE (*,*) 'SeedSize = ', SeedSize
      ALLOCATE(Seed(SeedSize))          ! 配列の割り当て
      DO I = 1, SeedSize
        CALL SYSTEM_CLOCK(count = Seed(i))! 現在時間を乱数の種に
      END DO
      CALL RANDOM_SEED(put = Seed(:))   ! 種を与える

      ! まずは求める真値を表示
      WRITE (*, '("Correct ",F13.7)') CORRECT
      ! モンテカルロ法開始
      I = 0
      K = 0
      DO WHILE(K .LT. Imax)
        DO I = 1, Istep
          ! (0,0) - (Ra, Rb) の中のランダムな点P
          CALL RANDOM_NUMBER(Px)
          CALL RANDOM_NUMBER(Py)
          Px = Px * Ra
          Py = Py * Rb

          ! 各焦点から、点P までの距離の合計
          Lab = SQRT((-F + Px) ** 2 + Py ** 2) + &
                SQRT(( F + Px) ** 2 + Py ** 2)
          IF (Lab .LE. L) THEN
              J = J + 1         ! 楕円の中ならカウントアップ
          END IF
          K = K + 1             ! 試行回数をカウント
        END DO
        ! 計算経過の表示
        CALL SYSTEM_CLOCK(Time_End)
        WRITE (*,100) RaRb4 * J / K, 100 * (1 - J/(K * AT_1)), K, &
                100.0*J/K, 100.0*K/Imax, DBLE(Time_End - Time_Start)/t_rate
      ENDDO

! 計算終了
! 真値 (再掲)
      WRITE (*, '("Correct ",F13.7)') CORRECT

      ! 終了時刻 実時間 & CPU 時間
      CALL SYSTEM_CLOCK(Time_End)
      CALL CPU_TIME(CpuTime_End)
      IF ( Time_End < Time_Start ) THEN
        Time_Diff = (t_max - Time_Start) + Time_End + 1
      ELSE
        Time_Diff = Time_End - Time_Start
      ENDIF

      Time_Diff = DBLE(Time_Diff/t_rate)        ! 単位を秒へ変換
      WRITE (*,'("elapsed time ", F0.3, "sec")') Time_Diff
      WRITE (*,'("elapsed cpu time ", F0.3, "sec")') &
        (CpuTime_End - CpuTime_Start)

! 面積 F13.7 (誤差、試行回数 I10、楕円の内側の割合、進行割合 F6.2%
100   FORMAT('  Aeria ', F13.7, ' diff:', F7.4, '% (Total:',I10 , &
        ' Inside:', F7.3,'% / ', F6.2'% Done / esp ', F0.2, ' sec)')

      END

いざ、実行！

tcsh

% gfortran EllipseArea_MonteCarlo_Solo.F90 && ./a.exe
 SeedSize =           33
Correct 45238.9354706
  Aeria 45227.4624000 diff: 0.0254% (Total:   1000000 Inside: 78.520% /   2.00% Done / esp 1.11 sec)
  Aeria 45235.6704000 diff: 0.0072% (Total:   2000000 Inside: 78.534% /   4.00% Done / esp 2.20 sec)
  Aeria 45239.8272000 diff:-0.0020% (Total:   3000000 Inside: 78.541% /   6.00% Done / esp 3.30 sec)
...snip...
  Aeria 45242.3640000 diff:-0.0076% (Total:  48000000 Inside: 78.546% /  96.00% Done / esp 52.45 sec)
  Aeria 45242.4125388 diff:-0.0077% (Total:  49000000 Inside: 78.546% /  98.00% Done / esp 53.55 sec)
  Aeria 45242.4142080 diff:-0.0077% (Total:  50000000 Inside: 78.546% / 100.00% Done / esp 54.64 sec)
Correct 45238.9354706
elapsed time 54.641sec (←時計の経過時間)
elapsed cpu time 54.547sec (← CPUを使っていた時間)

5千万回で 0.0077% って結構良いですねぇ。
まぁ、でもシングルスレッドだからねぇ。OpenMP にしたらもっと早いのかなぁ？