今更だけどモンテカルロ法で円周率推定

はじめに

先日これとほぼ同じ内容の記事を投稿しました。
しかしプログラムが不備だらけだったのでそちらは限定公開にします。
ぱっと思いついたものに書き換えたらうまくいったので投稿します。
もし気になる方がいれば　[問題の記事>>][3]
後で調べてみたのですが、今回のが割と正攻法みたいですね。

モンテカルロ法が分からない方は[この動画][2]で知ってください。

以上を踏まえて、

実行環境

Editor: VSCode
Shell: bash version 4.4.20
Compiler: gcc version 7.4.0

ソースコード

前回からの変更点は、
座標計算の方法と何となく処理時間が見れるようにしたこと。

Monte_Carlo2.c

# define _USE_MATH_DEFINES
# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>
# include <unistd.h>

# define POINT 940000       //１回の描画点数
# define LOOP 100.0         //モンテカルロ法実施回数

int main(){
    int in = 0, out = 0;
    double p = POINT, q = LOOP, ave = 0.0;
    for(int j = 0; j < q; j++) { //LOOP個のモンテカルロの円周率の平均値を求める
        int sum = 0;
        srand(time(NULL));
        for(int i = 0; i < p; i++){
            double x = (double)rand() / RAND_MAX;
            double y = (double)rand() / RAND_MAX;
            if(sqrt(pow(x, 2) + pow(y, 2)) <= 1){
                //printf("内[%1.4f, %1.4f]\n", x, y);
                /*↑コメントアウトを開放すると大量に座標見れる*/
                sum++;
                in++;       //四分円内にあった回数
            }
            else {
                //printf("　[%1.4f, %1.4f]外\n", x, y);
                /*↑コメントアウトを開放すると大量に座標見れる*/
                out++;      //四分円外にあった回数
            }
        }
        ave += sum / p;
    }
    double run = (double)clock() / CLOCKS_PER_SEC;
    double pi = 4 * ave / q;
    printf("処理時間：%f[秒]\n", run);
    printf("１回の描画点数：%d[個]\n", (int)p);
    printf("モンテカルロ法：%d[回]\n", (int)q);
    printf("四分円内の点数：%d[個]\n", in);
    printf("四分円外の点数：%d[個]\n", out);
    printf("円周率：%1.7f\n", pi);
    printf("今回の精度：%2.5f％\n", (M_PI-fabs(M_PI-pi))/M_PI*100);

    return 0;
}

実行方法, 実行結果

-lmは忘れずにつけましょう。

bash

$ gcc -Wall Monte_Carlo2.c -lm
$ ./a.out
処理時間：6.703125[秒]
１回の描画点数：940000[個]
モンテカルロ法：100[回]
四分円内の点数：73826428[個]
四分円外の点数：20173572[個]
円周率：3.1415501
今回の精度：99.99865％

おわりに

描画点数、モンテカルロ法の実施回数色々変えて暇つぶしてください。
よいプログラミングライフを(送りたい)。

参考

• 一瞬ルートでつまずいた：[ライブラリのリンクを忘れずに][1]
• 作るきっかけになった動画：[【すげえw】モンテカルロ法による中学生も理解できる円周率の求め方が面白い【物理エンジン】][2]
• 不備だらけの前投稿：[今更だけどモンテカルロ法で円周率推定][3]

スペシャルサンクス：@fujitanozomu
[1]:https://qiita.com/percipere/items/678a8ec1299a7564f6a3
[2]:https://youtu.be/dn-RtqY9aCQ
[3]:https://qiita.com/RoaaaA/private/38ef42576985f67be720