電子ブロック工房 IC・トランジスタで出来たコンピューターを設計・製作するためのブログ

モンテカルロ法で円周率を求めるアルゴリズムはいたって簡単。-1 から 1 までの乱数を、x, y として用意し、図のように点(x,y)と座標の原点までの距離が１以下であれば円の内側だとみなされる。ピタゴラスの定理より、単に


の値が１より大きいか小さいかを求めればよい。これを何回も繰り返せば、円の内側に点が入る確率が求まる。その確率を p とすると、外側の正方形の面積が４であるから、円周率は次のように得られる。



これを、C でシミュレーションするには、以下の様な簡単なコードでよい。

#include <stdlib.h>
#include <stdio.h>
#include <conio.h>

float sub();

void main (){
    for(unsigned int i=0;!kbhit();i++){
        float pai=sub();
        printf( "%d %f\n", i, pai);
    }
}

float sub(){
    static float num,inside;
    float x,y,pai;
    float max=RAND_MAX+1;
    
    // Decide the position randomly
    x=((float)rand()+(float)rand()/(float)max)/(float)max;
    y=((float)rand()+(float)rand()/(float)max)/(float)max;
    
    // Determine if inside the circle
    num++;
    if (x*x+y*y<1) inside++;

    // Calculate the pai
    pai=4*inside/num;
    
    return pai;
}

プログラム簡略化のため、0≦x≦1, 0≦y≦1 としている。１０万回ほど計算すれば、3.14 が求められる。ちなみに、このアルゴリズムでは 3.141 まで求められたが、3.1416 は出せなかった。float の許容範囲を超えているようである。

浮動小数点を、double で使用すると…

#include <stdlib.h>
#include <stdio.h>
#include <conio.h>

double sub();

void main (){
    for(unsigned int i=0;!kbhit();i+=10000){
        double pai=sub();
        printf( "%d %f\n", i, pai);
    }
}

double sub(){
    static double num,inside;
    double x,y,pai;
    double max=RAND_MAX+1;
    
    for(int i=0;i<10000;i++){
        // Decide the position randomly
        x=((double)rand()+(double)rand()/(double)max)/(double)max;
        y=((double)rand()+(double)rand()/(double)max)/(double)max;
        
        // Determine if inside the circle
        num++;
        if (x*x+y*y<1) inside++;
    }

    // Calculate the pai
    pai=4*inside/num;
    
    return pai;
}

この場合は、５億回ほど計算すると、3.1416 まで求めることが出来た。この辺りが、一つの目標であろう。

　Ｃでさらに計算を続けると、４０億回の計算で、3.14159 まで出せた。これより下の桁を出すには、乱数の精度を上げる必要があるかもしれない。