例 : 円周率の計算

ここまで学んだことを基に例題を見てみる。ここでは疑似乱数を用いて円周率を求めるルーチンを作成する。

今、半径 1 の円と、それに外接する正方形を考える。それぞれの面積は

$$\begin{eqnarray*} S_{正方形} = (2\times 1)^2 = 4 \\ S_{円} = \pi \times 1^2 = \pi \end{eqnarray*}$$

である。さて、この正方形の中にランダムに点を打っていくと、正方形内に打ち 込まれた点の数 $N_{正方形}$ と円内に打ち込まれた点の数 $N_{円}$ の比はそれぞれの面積の比に等しくなる。

$$\begin{eqnarray*} N_{円}/N_{正方形} = S_{円}/S{正方形} \end{eqnarray*}$$

これと先の式を合わせると、

$$\begin{eqnarray*} \pi = 4N_{円}/N_{正方形} \end{eqnarray*}$$

となり、$\pi$ の値を求めることができる。

さて、プログラムを見てみよう。²

ludolph.c ------------------------------------------------------------

#include<stdio.h>

int main(){

  int NEVENT = 10000;                    /* 正方形に入る点の数 */
  int IA = 421,IC = 54773,IM = 259200;   /* 乱数のパラメタ */
  int iseed = 12345;                     /* 乱数のたね */
  int n = 0;                             /* 円内に入った点の数。初期値は 0 */
  int i,j;                               /* for 文のカウンタ */
  double x,r, pi;                        /* x(y)座標, 動径の大きさ, 円周率 */

  for(i = 0; i < NEVENT; i++){           /* NEVENT だけループを回す */
    r = 0.0;
    for(j = 0; j < 2; j++){ 
                    /* 2 回まわすことで、変数 x に x、y 両方の役割を負わせる */
      iseed = (iseed*IA+IC)%IM;          /* 乱数発生ルーチン */
      x = (double)iseed/(double)IM;
      r += x*x;                          /* r = r+x*x と同じ */
    }
    if(r <= 1.0){               /* r が1より小さいか同じなら n に 1 を加える */
      n++;
    }
  }
  pi = 4.0*(double)n/(double)NEVENT;     /* n : NEVENT = pi:4.0 */
  printf("pi = %f \n",pi);               /* 結果の表示 */
}

----------------------------------------------------------------------

上から順に見ていく。

メイン関数全体を表している。メイン関数の上にある、

#include<stdio.h>

については後で説明する。

変数の定義を行っている。変数のプログラム内での意味は後で説明する。ここでは先程述べなかった点が2つ含まれている。

• 初期値の代入 : 変数に初期値を代入したい時は、変数名と代入する値を ''='' で結ぶ。
• 変数名の規則 : と言うと大袈裟だが、大文字と小文字は区別される。''a'' と ''A'' は別の変数と見なされる。通常、全て大文字で書いた場合は定数とする(#define 参照)。

for 文を用いて正方形内に入る点の数だけループをまわす。つまり、変数 i が 0 〜 NEVENT-1 までの計 NEVENT 回中の実行文を行う。この中で、点が円の中に入ったか判定し、入ったらその数を数えるという行為を繰り返す。

ここでは $(x,y)$ の組を決めるのだが、ここでは $x$ と $y$ をそれぞれ別の変数に代入するのではなく、一つの変数 x に2つの値の役割を負わせ、その自乗和を計算することにすることにした。

まずこの部分から見ていく。ここで新たに出てきた記号 ''+='' は、''r = r+x*x'' 、すなわち「r+x*x を計算し、それを再び r に代入する」ことを意味する。このような使い方は四則演算の記号全てに可能である。

a += b -> a = a + b

a -= b -> a = a - b

a *= b -> a = a * b

a /= b -> a = a / b

a %= b -> a = a % b

そうしてみると、まず r を初期化してから for 文に入る。 for 文内では、まず変数 x に座標 $x$ の役を負わせ、 r に x の自乗を足す。次に再び前に戻って、変数 x に座標 $y$ の役を負わせて、r に x の自乗を 足す。2 回まわったので、for 文の条件から外れて、for 文の外へ出る。こうして動径方向の自乗和が求まる。

この2行で疑似乱数を生成する。乱数発生法はここでは省略する³。ここで注目して欲しいのは ''x = (double)iseed/(double)IM;'' という計算法である。変数 iseed, IM は int 型に対し、変数 x は double 型である。この場合、通常通り ''x = iseed/IM;'' と計算すると、iseed と IM が int 型のため、iseed/IM を int 型の答えで一旦出した後に double 型に変換して x に代入される。そのため、x には 0 (まれに 1 )しか入らない。そのためここではこれら int 型の変数の前に ''(double)'' と付けて double 型とみなして計算し、計算結果も double 型で求める。そうすることで、x には 0 〜 1の間の任意の実数が代入される。

この部分は 「3.4」 で挙げた例と同じである。

ここで上の関係式を用いて $\pi$ を導出する。

ここでは結果を画面に出力するために、標準関数 printf を用いている。標準関数とは OS であらかじめ用意している関数のことである。標準関数は種類ごとにパッケージにされていて(標準ライブラリと言う)、使用する標準関数に合わせてメイン関数の前に指定する必要がある。

printf の場合、標準ライブラリは stdio.h を指定すれば良い。指定の仕方は、ソースの一番はじめに

#include<stdio.h>

とおけば良い。

printf の一般的な書式は

  printf("データの書式や出力させる文字列等  … \n", 出力される変数, …)

で表される。変数内の値を出力したい時は下の表のような記号をいれる。'\n' は改行を表す。それ以外は、そのまま出力される。

%f	実数
%d	整数
%s	文字列

この例では、画面には

pi = (値)

の様に出力される。