Jeans & Development コンピューターのことなどを綴ったメモ (旧：目から鱗 w/SQLite)

ちょっと訳有って、円周率（π、パイ、pi）を、コンピューターに計算させてみた。とりあえず、例によって、PHPで。

円周率の計算は、手作りコンピューターでやろうとしていることなので、アルゴリズムとしては最も簡単なものを利用した。最初に作成したスクリプトは、以下のとおり。

for ($max=1;$max<=8192;$max=$max*2) {
    $t=time();
    echo "$max: ".calc_pi($max)."  (".(time()-$t).")\n";
}

function calc_pi($max) {
    $total=0;
    for ($x=0;$x<$max;$x++) {
        for ($y=0;$y<$max;$y++) {
            if ($x*$x + $y*$y < $max * $max) $total++;
        }
    }
    $pi=4 * ( $total/($max*$max) );
    return $pi;
}

正方形の内部にまんべんなく点を打ったとき、その点から一つの頂点までの距離が正方形の辺の長さより短いかどうかをピタゴラスの定理を用いて判断することで、円の内部かどうかを調べ、円の面積から円周率を計算している。

実行結果は、以下のようになった。

1: 4  (0)
2: 4  (0)
4: 3.75  (0)
8: 3.5  (0)
16: 3.34375  (0)
32: 3.25390625  (0)
64: 3.19921875  (0)
128: 3.170166015625  (0)
256: 3.15673828125  (0)
512: 3.1490936279297  (0)
1024: 3.1453971862793  (1)
2048: 3.1435289382935  (2)
4096: 3.1425535678864  (10)
8192: 3.1420764327049  (40)

上のスクリプトは正方形内のすべての点を調べており、あまりにも効率が悪い。少し修正して、以下のようなスクリプトにした。

for ($max=1;$max<=536870912;$max=$max*2) {
    $t=time();
    echo "$max: ".calc_pi($max)."  (".(time()-$t).")\n";
}

function calc_pi($max) {
    $total=0;
    $max2=$max*$max;
    $prev_y=$max-1;
    for ($x=0;$x<$max;$x++) {
        $max2_x2=$max2-$x*$x;
        for ($y=$prev_y;0<=$y;$y--) {
            if ($y*$y < $max2_x2) break;
        }
        $total+=$y+1;
        $prev_y=$y;
    }
    $pi=4 * ( $total/$max2 );
    return $pi;
}

実行結果は、以下のとおり。

1: 4  (0)
2: 4  (0)
4: 3.75  (0)
8: 3.5  (0)
16: 3.34375  (0)
32: 3.25390625  (0)
64: 3.19921875  (0)
128: 3.170166015625  (0)
256: 3.15673828125  (0)
512: 3.1490936279297  (0)
1024: 3.1453971862793  (0)
2048: 3.1435289382935  (0)
4096: 3.1425535678864  (0)
8192: 3.1420764327049  (0)
16384: 3.1418349444866  (0)
32768: 3.1417138241231  (0)
65536: 3.1416533300653  (0)
131072: 3.1416230569594  (0)
262144: 3.1416078822222  (1)
524288: 3.1416002743645  (0)
1048576: 3.1415964659827  (2)
2097152: 3.141594559359  (2)
4194304: 3.1415936068045  (6)
8388608: 3.1415931302873  (12)
16777216: 3.141592891993  (22)
33554432: 3.1415927727818  (46)
67108864: 3.1415927131891  (91)
134217728: 3.1415926857583  (184)
268435456: 3.141592670158  (365)
536870912: 3.1415926618737  (730)

小数点以下７桁目を確定させるのに、約３分。まあ、こんなものかな。後で、別のプログラミング言語でも調べてみることにしよう。