モンテカルロ法の面積近似

Mが2の32乗の、線形合同法で乱数を発生し、
モンテカルロ法を使って、
領域a = {(x,y)|x≧1,y≧2,(x-1)^2 + (y-2)^2 ≦1}の面積の近似値を計算するプログラムを作成し、これを利用してn=10000とした場合の上の領域の面積、および円周率πの近似値を求めよ。
(関数を使う)

という課題が出たのですが、乱数の出し方？が変であまり近似されません(^_^;)
どなたかどこが変か教えてください(>_<)

これだと、2と4になってしまうんです(/_;)

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

#define N 100

/*領域Rを含む面積T　領域R=求める面積*/
#define Nx 1
#define Ny 2

/*乱数を発生させる関数*/
unsigned long int random_g(unsigned long int x0, unsigned long int a, unsigned long int c);

int main(){
unsigned long int x0, a, c, M;

double x, y, s;
int n, count_in, count_out;

count_in = 0;
count_out = 0;

printf("x0: ", x0 ); scanf("%d", &x0);

a = 61;
c = 49;

for(n = 0; n < 10000; n++){
/*乱数を発生*/

x = (random_g(x0, 4*a*n+1, c) / (double)4294967296)*Nx + 1;
y = (random_g(x0, 4*a*n+1, c) / (double)4294967296)*Ny + 1;

if(x >= 1 && y >= 1 && pow(x-1, 2) + pow(y-2, 2) <= 1)
count_in++;
else
count_out++;
}

s = (double)count_in/(count_in + count_out)*Nx*Ny;
printf("s = %f\n", s);
printf("pi = %f\n", 2*s);


return 0;
}

unsigned long int random_g(unsigned long int x0, unsigned long int a, unsigned long int c){

int i;

x0 = (a*x0 + c);
/*printf("%u\n", x0);*/

return x0;
}

No.3 ベストアンサー 回答者： asuncion 回答日時： 2007/04/30 00:43

正しいかどうかよくわからないのですが、こんな感じなのでしょうか。

求める面積は(1, 2)を中心とする半径1の円の4分の1ですから、
円周率の値は面積を4倍しなければなりません。

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

#define N (10000)
#define Nx (1)
#define Ny (2)

/*乱数を発生させる関数*/
unsigned long random_g(unsigned long x0, unsigned long a, unsigned long c);

int main(void)
{
　　unsigned long x0, a, c;
　　double x, y, s;
　　int n, count_in;
　　
　　printf("x0: ");
　　scanf("%d", &x0);
　　a = 61;
　　c = 49;
　　
　　for (count_in = n = 0; n < N; n++) {
　　　　/*乱数を発生*/
　　　　x0 = random_g(x0, a, c);
　　　　x = x0 / (double) 4294967296 + Nx;
　　　　x0 = random_g(x0, a, c);
　　　　y = x0 / (double) 4294967296 + Ny;
　　　　if (x >= Nx && y >= Ny && (x-Nx) * (x-Nx) + (y-Ny) * (y-Ny) <= 1)
　　　　　　count_in++;
　　}
　　
　　s = (double) count_in / n;
　　printf("s = %f\n", s);
　　printf("pi = %f\n", 4 * s);
　　return 0;
}

unsigned long random_g(unsigned long x0, unsigned long a, unsigned long c)
{
　　return a * x0 + c;
}

（注）インデントのため、全角空白を使っています。

この回答へのお礼

ありがとうございます！
参考になりました！！
でも、
/*乱数を発生させる関数*/
unsigned long random_g
としているので、乱数はrandom_g内で発生させたいのですが、
それは無理ですか？？(T_T)

お礼日時：2007/04/30 00:53

No.4 回答者： asuncion 回答日時： 2007/04/30 09:06

> unsigned long random_g

> としているので、乱数はrandom_g内で発生させたい

random_g()の中で、前回発生させた乱数(x0)をもとに、
所定のパラメータ(a, c)を用いて新しい乱数(returnする値)を
発生させています。

この回答へのお礼

そうですねっっ！！
勘違いしてたみたいです。。

ありがとうございました(^_^)

お礼日時：2007/04/30 19:18

No.2 回答者： asuncion 回答日時： 2007/04/29 22:14

> 領域a = {(x,y)|x≧1,y≧2,(x-1)^2 + (y-2)^2 ≦1}

この定義と
> if(x >= 1 && y >= 1 && pow(x-1, 2) + pow(y-2, 2) <= 1)
このif文とが一致していませんが、よいのでしょうか？

また、このif文の真偽によってcount_inかcount_outの
いずれか一方だけをインクリメントするのですから、
count_inとcount_outの和は試行回数nと一致しますよね。
よって、count_outは特に必要ないと思います。

この回答への補足

確かにそうですね(^_^;)

補足日時：2007/04/30 00:49

No.1 回答者： rabbit_cat 回答日時： 2007/04/29 21:51

まずは、random_gの中でせっかく計算した、x0を覚えておくようにしないといけません。

そうですね。

１　x0をグローバル変数にする。
２
　x = (random_g(x0, 4*a*n+1, c) / (double)4294967296)*Nx + 1;
　を
　x = (random_g(x0, a, c) / (double)4294967296)*Nx + 1;
　に変更する。
　yについても同様に、4*a*n+1をaに変更する。

とどうでしょうかね。
２で指摘した4*a*n+1ていうのは、apple_cubeさんが何を考えてそうなっていたのか意図をつかみかねているので、もしかしたら、方向違いの回答になっているかもしれません。

この回答へのお礼

ありがとうございます！！
4*a*n+1は乱数の数が近くてどうすれば全然関係なくなるのか
考えて色々試していた残骸です…(^_^;)

それで、x0をグローバルか変数にしたのですが、答えがやはり2と4になってしまいました(>_<)

お礼日時：2007/04/29 22:04