* 円周率の計算を最適化してみよう！
-(by [[K]], 2015.11.10)
*** はじめに
-ここでは原則としてプログラムの話題は扱わない方針だったのだけど、[[math/03]]の内容に関連していることもあり、これはあえてここで扱うことにする。
*** 改良前の確認
-以下が改良前のプログラムになる。
#include <stdio.h>
#define N 10000
int main()
{
int i, j, k, c[5];
double x0, y0, x1, y1;
for (k = 0; k < 5; k++)
c[k] = 0;
for (j = 0; j < N; j++) {
y0 = (double) j / N;
y1 = (double) (j + 1) / N;
for (i = 0; i < N; i++) {
x0 = (double) i / N;
x1 = (double) (i + 1) / N;
k = 0;
if (x0 * x0 + y0 * y0 <= 1.0) k++;
if (x1 * x1 + y0 * y0 <= 1.0) k++;
if (x0 * x0 + y1 * y1 <= 1.0) k++;
if (x1 * x1 + y1 * y1 <= 1.0) k++;
c[k]++;
}
}
double pi0 = 4.0 * c[4] / ((double) N*N);
double pi1 = 4.0 * (c[1] + c[2] + c[3] + c[4]) / ((double) N*N);
printf("pi=%.8f\n", (pi0 + pi1) / 2);
return 0;
}
-このプログラムは[[math/03]]のプログラムを少しだけ改良したもので、やっていることはパイの上限下限を推定して、その単純平均を出力しているだけだ。
-このプログラムは、少しがんばるとみるみる小さくできるので、それを見せたいと思う。
***左下と右上しかチェックしない
-プログラムでは、c[4]の値とc[1〜4]の合計値しか必要としていない。これを判定するだけならif文は2つで十分だ。つまり右上が円盤の中に入っていればミニ正方形が1個とカウントされるべきだし、左下だけが円盤の中に入っていればミニ正方形が0.5個とカウントされるべきだ。
-0.5とかだとめんどくさいので、2倍で数えることにして、最後の計算では4倍ではなく2倍にすることにしよう。
#include <stdio.h>
#define N 10000
int main()
{
int i, j, c = 0;
double x0, y0, x1, y1;
for (j = 0; j < N; j++) {
y0 = (double) j / N;
y1 = (double) (j + 1) / N;
for (i = 0; i < N; i++) {
x0 = (double) i / N;
x1 = (double) (i + 1) / N;
if (x0 * x0 + y0 * y0 <= 1.0) c++;
if (x1 * x1 + y1 * y1 <= 1.0) c++;
}
}
double pi = 2.0 * c / ((double) N*N);
printf("pi=%.8f\n", pi);
return 0;
}
***if文を一つにする
-i=0やj=0のような辺のケースを除いて、任意の格子点に注目すると、それはミニ正方形の左下の場合と右上の場合とで、2度評価されている。これは無駄だ。一つの格子点が一度だけ評価されるほうが高速である。
-これだと++が2回実行されていたのが1回に減ってしまうので、値を2倍してやる。
-そして辺の部分の寄与については、2N-1を足してやればよい。
-これで値が全く変わることなく、速度は2倍くらいになる。プログラムも短くなる。
#include <stdio.h>
#define N 10000
int main()
{
int i, j, c = 0;
double x, y;
for (j = 1; j < N; j++) {
y = (double) j / N;
for (i = 1; i < N; i++) {
x = (double) i / N;
if (x * x + y * y <= 1.0) c++;
}
}
c *= 2;
c += 2 * N - 1;
double pi = 2.0 * c / ((double) N*N);
printf("pi=%.8f\n", pi);
return 0;
}
***整数演算だけにする
-もし最後の計算結果が(N*N)倍になってもよければ、doubleなんか使わなくても計算できる。今N=1万を考えているので、N*Nはちょうど1億になる。結果が1億倍になっても、パイの値が3.14くらいであるということはみんな知っているので、どこに小数点を打てばいいのか混乱することはまずない。ということで、doubleを使わずに整数演算だけにしてしまおう。
#include <stdio.h>
#define N 10000
int main()
{
int i, j, c = 0;
for (j = 1; j < N; j++) {
for (i = 1; i < N; i++) {
if (i * i + j * j <= N * N) c++;
}
}
printf("pi=%d\n", c * 4 + 4 * N - 2);
return 0;
}
-おお、なんということか！ こんな単純なプログラムでも円周率は計算できてしまうのだ。
***高速化
-上記のプログラムは、当然ながらif文が(N*N)回実行される。これは遅い。・・・それならば高速化しようではないか。
-上記のプログラムは、当然ながらif文が(N*N)回実行される。これはまだ遅い。・・・それならば高速化しようではないか。
#include <stdio.h>
#define N 10000
int main()
{
int i = N - 1, j, c = 0;
for (j = 1; j < N; j++) {
while (i * i + j * j > N * N) i--;
c += i;
}
printf("pi=%d\n", c * 4 + 4 * N - 2);
return 0;
}
-これでwhileの中の条件判定は(N*2)回くらいしか評価されないことになる。やっていることは、cを何度も++するのではなくて、横一列をどかっとまとめて足しているだけである。
***まとめ
-いきなり最後のプログラムを示されても理解は難しいと思うが、このように順を追えば、そこそこ分かってもらえるのではないかと思う。
-整数演算だけで、しかも割り算や無理関数を全く使わずとも、この程度のプログラムで円周率は計算できるのだ。
*こめんと欄
#comment
![[K Wiki]](k.jpg "[K Wiki]")
