利用蒙特卡洛算法求圆周率是一个概率的方法,关于这方面的内容很多,而且也很容易理解,更多具体分析过程可以参考如下文章:
http://blog.csdn.net/zzhflyqh/article/details/2288722
http://www.oschina.net/code/snippet_76_4482
http://yangyi.iteye.com/blog/996177
下面是我的理解和代码:
蒙特卡洛算法是通过概率来计算pi的值的。对于一个单位为1的正方形,以其某一个顶点为圆心,边为半径在正方形内画扇形(一个1/4的圆形的扇形),那么扇形的面积就是pi/4。这样,利用概率的方式,“随机”往正方形里面放入一些“点”,根据这些点在扇形内的概率(在扇形内的点数/投的总点数),就可以得到扇形的面积。
简单理解,由于正方形为单位1的边长,所以正方形面积S1=1. 假设随机投如N个点,其中M个落入扇形内,那么可以认为扇形的面积为M/N*1。当然,本身,这就是一个概率的过程,不是精确的。当N足够大的时候,得到的面积更精确。由于已经扇形的面积为pi/4,所以根据M/N*1=pi/4就可以计算出pi的值。
这里的关键问题是:(1)如何随机投入点?由于C语言提供了随机函数,所以利用随机函数就可以模拟随机投点的过程。(2)如何判断点在扇形内?这时候,我们先考虑建立的坐标系,假设单位正方形的顶点在原点,正方形在第一象限,即正方形范围为(0,0)到(1,1)的坐标内,并且假设我们扇形的圆心也为原点,那么根据圆的方程,假设投的点坐标为(x,y),如果x*x+y*y<1,则点(x,y)在扇形内。
参考下图进行理解:
pi = 4*M/N
下面是具体的程序:
#include #include #include #include #define N 10000000000
#define int_t __int64 // the type of N and M
// Rand from 0.0 to 1.0
inline double rand01()
{
return rand()*1.0/RAND_MAX;
}
int main()
{
srand((unsigned)time(NULL));
double x,y;
int_t M = 0;
for (int_t i = 0;i < N;i++)
{
x = rand01();
y = rand01();
if (x*x+y*y<1) M++;
}
double pi = (double)4*M/N;
printf("%lf\n", pi);
}
写这个例子的目的是这是一个适合用多线程处理的很好的例子,在进行其它的多线程实现的时候,将不再分析这个过程。