#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
 
#define NUM_POINTS 1000000 // 最大乱数の個数
#define TRUE_AREA (M_PI - 3 * sqrt(3) + 3) / 3 // 真の面積
#define NUM_STEPS 12 // 横軸の乱数の個数を12個に制限
 
// 2つの円の内部にあるかどうかを判定する関数
int is_inside(double x, double y) {
    double r = 1.0;
    // 円1：中心(0,1)、半径1
    if ((x - 0) * (x - 0) + (y - 1) * (y - 1) <= r * r && 
        // 円2：中心(1,0)、半径1
        (x - 1) * (x - 1) + (y - 0) * (y - 0) <= r * r) {
        return 1;
    }
    return 0;
}
 
int main() {
    int i, inside_count = 0;
    double x, y;
    double estimated_area;
    double error;
    int step_size = NUM_POINTS / NUM_STEPS; // 乱数の個数を12ステップに分割
 
    // 乱数の生成を初期化
    srand(time(NULL));
 
    // タイトルを出力（列名として表示される）
    printf("乱数の個数,誤差\n");
 
    // Monte Carlo法で面積を近似する
    for (i = 1; i <= NUM_POINTS; i++) {
        // 0から1の範囲で乱数を生成
        x = (double)rand() / RAND_MAX;
        y = (double)rand() / RAND_MAX;
 
        // 円の内部にあるかどうかを判定
        if (is_inside(x, y)) {
            inside_count++;
        }
 
        // 12等分されたステップごとに誤差を計算して出力
        if (i % step_size == 0) {
            estimated_area = 4.0 * (double)inside_count / i; // 4倍するのは単位正方形全体の面積が1×1なので
            error = fabs(estimated_area - TRUE_AREA);
 
            // カンマ区切りで乱数の個数と誤差を出力
            printf("%d,%f\n", i, error);
        }
    }
 
    printf("計算が完了しました。\n");
 
    return 0;
}
 

I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxzdGRsaWIuaD4KI2luY2x1ZGUgPG1hdGguaD4KI2luY2x1ZGUgPHRpbWUuaD4KCiNkZWZpbmUgTlVNX1BPSU5UUyAxMDAwMDAwIC8vIOacgOWkp+S5seaVsOOBruWAi+aVsAojZGVmaW5lIFRSVUVfQVJFQSAoTV9QSSAtIDMgKiBzcXJ0KDMpICsgMykgLyAzIC8vIOecn+OBrumdouepjQojZGVmaW5lIE5VTV9TVEVQUyAxMiAvLyDmqKrou7jjga7kubHmlbDjga7lgIvmlbDjgpIxMuWAi+OBq+WItumZkAoKLy8gMuOBpOOBruWGhuOBruWGhemDqOOBq+OBguOCi+OBi+OBqeOBhuOBi+OCkuWIpOWumuOBmeOCi+mWouaVsAppbnQgaXNfaW5zaWRlKGRvdWJsZSB4LCBkb3VibGUgeSkgewogICAgZG91YmxlIHIgPSAxLjA7CiAgICAvLyDlhoYx77ya5Lit5b+DKDAsMSnjgIHljYrlvoQxCiAgICBpZiAoKHggLSAwKSAqICh4IC0gMCkgKyAoeSAtIDEpICogKHkgLSAxKSA8PSByICogciAmJiAKICAgICAgICAvLyDlhoYy77ya5Lit5b+DKDEsMCnjgIHljYrlvoQxCiAgICAgICAgKHggLSAxKSAqICh4IC0gMSkgKyAoeSAtIDApICogKHkgLSAwKSA8PSByICogcikgewogICAgICAgIHJldHVybiAxOwogICAgfQogICAgcmV0dXJuIDA7Cn0KCmludCBtYWluKCkgewogICAgaW50IGksIGluc2lkZV9jb3VudCA9IDA7CiAgICBkb3VibGUgeCwgeTsKICAgIGRvdWJsZSBlc3RpbWF0ZWRfYXJlYTsKICAgIGRvdWJsZSBlcnJvcjsKICAgIGludCBzdGVwX3NpemUgPSBOVU1fUE9JTlRTIC8gTlVNX1NURVBTOyAvLyDkubHmlbDjga7lgIvmlbDjgpIxMuOCueODhuODg+ODl+OBq+WIhuWJsgoKICAgIC8vIOS5seaVsOOBrueUn+aIkOOCkuWIneacn+WMlgogICAgc3JhbmQodGltZShOVUxMKSk7CgogICAgLy8g44K/44Kk44OI44Or44KS5Ye65Yqb77yI5YiX5ZCN44Go44GX44Gm6KGo56S644GV44KM44KL77yJCiAgICBwcmludGYoIuS5seaVsOOBruWAi+aVsCzoqqTlt65cbiIpOwoKICAgIC8vIE1vbnRlIENhcmxv5rOV44Gn6Z2i56mN44KS6L+R5Ly844GZ44KLCiAgICBmb3IgKGkgPSAxOyBpIDw9IE5VTV9QT0lOVFM7IGkrKykgewogICAgICAgIC8vIDDjgYvjgokx44Gu56+E5Zuy44Gn5Lmx5pWw44KS55Sf5oiQCiAgICAgICAgeCA9IChkb3VibGUpcmFuZCgpIC8gUkFORF9NQVg7CiAgICAgICAgeSA9IChkb3VibGUpcmFuZCgpIC8gUkFORF9NQVg7CgogICAgICAgIC8vIOWGhuOBruWGhemDqOOBq+OBguOCi+OBi+OBqeOBhuOBi+OCkuWIpOWumgogICAgICAgIGlmIChpc19pbnNpZGUoeCwgeSkpIHsKICAgICAgICAgICAgaW5zaWRlX2NvdW50Kys7CiAgICAgICAgfQoKICAgICAgICAvLyAxMuetieWIhuOBleOCjOOBn+OCueODhuODg+ODl+OBlOOBqOOBq+iqpOW3ruOCkuioiOeul+OBl+OBpuWHuuWKmwogICAgICAgIGlmIChpICUgc3RlcF9zaXplID09IDApIHsKICAgICAgICAgICAgZXN0aW1hdGVkX2FyZWEgPSA0LjAgKiAoZG91YmxlKWluc2lkZV9jb3VudCAvIGk7IC8vIDTlgI3jgZnjgovjga7jga/ljZjkvY3mraPmlrnlvaLlhajkvZPjga7pnaLnqY3jgYwxw5cx44Gq44Gu44GnCiAgICAgICAgICAgIGVycm9yID0gZmFicyhlc3RpbWF0ZWRfYXJlYSAtIFRSVUVfQVJFQSk7CgogICAgICAgICAgICAvLyDjgqvjg7Pjg57ljLrliIfjgorjgafkubHmlbDjga7lgIvmlbDjgajoqqTlt67jgpLlh7rlipsKICAgICAgICAgICAgcHJpbnRmKCIlZCwlZlxuIiwgaSwgZXJyb3IpOwogICAgICAgIH0KICAgIH0KCiAgICBwcmludGYoIuioiOeul+OBjOWujOS6huOBl+OBvuOBl+OBn+OAglxuIik7CgogICAgcmV0dXJuIDA7Cn0K