蒙特卡罗法求蒲丰(buffon)问题

// 蒙特卡罗法_蒲丰氏问题.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include "stdlib.h" #include "time.h" #include "math.h" #define g_nNumOfExperiment 10000 #define TINY 0.000001 #define MAX_ITER 10000 //本程序用蒙特卡罗法求解蒲丰氏问题 int main(int argc, char* argv[]) { int i; double x[g_nNumOfExperiment+1];//针的中点所在的x轴坐标 double ang[g_nNumOfExperiment+1];//针与平行直线的夹角 double total; double f_a=4.0;//平行线的间距的一半 double f_l=3.0;//针的长度的一半 double pi=0.; double pa;//针与平行线相交的频率 double sum=0.; printf("蒙特卡罗法求解蒲丰问题\n"); for(int iter=1;iter<=MAX_ITER;iter++) { total=0.; /* Seed the random-number generator with current time so that * the numbers will be different every time we run. */ srand( (unsigned)time( NULL ) ); /* Generate random numbers. */ for( i = 1; i <= g_nNumOfExperiment;i++ ) { x[i]=double(rand()); x[i]=x[i]*f_a/RAND_MAX; ang[i]=double(rand()); ang[i]=f_l*sin(ang[i]);//这里认为夹角在0到2*pi中均匀分布 //最后将结果除以2,即可得到pi if (ang[i]>=x[i]) { ang[i]=1.; } else ang[i]=0.; total+=ang[i]; } pa=total/g_nNumOfExperiment;//针与平行线相交的频率 pi=2*f_l/(f_a*pa);//求圆周率 sum+=pi; } pi=sum/MAX_ITER; //输出结果 printf("圆周率为:%f\n",pi/2); return 0; }