【多核程序设计】蒙特卡洛算法求pi

背景知识：

蒙特卡洛是摩纳哥公国第一大城市，与澳门、美国拉斯维加斯并称世界三大赌城。位于地中海沿岸，首都摩纳哥之北，建于阿尔卑斯山脉突出地中海的悬崖之上。景色优美，是地中海地区旅游胜地。市内建有豪华的旅馆、俱乐部、歌剧院、商店、游泳池、温泉浴室、运动场等娱乐设施 。城内开设有蒙特卡洛大赌场。赌场建于1865年，为双层楼建筑，上有钟楼、塔厅和拱形亭阁，还饰以若干人物雕塑，庭前棕榈树成行，还辟有花园，旁边有大酒店和酒吧间。整个城市在旺季时，约有赌场70多个，约有赌室3500间左右。蒙特卡罗赌场由国家经营。当地的其他活动，许多也带有赌博色彩。游客住的旅店房间，有抽奖的号码，中奖的免付部分房费。早餐的牛奶麦片粥里，如遇上金属牌子 ，亦可领奖。该城只有1万人口，但每天报纸销量可达100万份 ，因为报纸上都印有可能得奖的号码。游客最后离境，购买的车票上也印有彩票号码，于离境前开彩。经营赌业是摩纳哥的主要经济来源，每年都从赌业中收取高额外汇利润。

蒙特卡洛算法简单描述：

以概率和统计理论方法为基础的一种计算方法。将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解。比如，给定x=a，和x=b，你要求某一曲线f和这两竖线，及x轴围成的面积，你可以起定y轴一横线 y=c 其中c>=f(a) and c>=f(b)，很简单的，你可以求出 y=c,x=a,x=b,及 x轴围成的矩形面积，然后利用随机参生生大量在这个矩形范围之类的点，统计出现在曲线上部点数和出现在曲线下部点的数目，记为：doteUpCount,nodeDownCount,然后所要求的面积可以近似为 doteDownCounts所占比例*矩形面积。

问题描述：

在数值积分法中，利用求单位圆的1/4的面积来求得Pi/4从而得到Pi。单位圆的1/4面积是一个扇形，它是边长为1单位正方形的一部分。只要能求出扇形面积S1在正方形面积S中占的比例K=S1/S就立即能得到S1，从而得到Pi的值。怎样求出扇形面积在正方形面积中占的比例K呢？一个办法是在正方形中随机投入很多点，使所投的点落在正方形中每一个位置的机会相等看其中有多少个点落在扇形内。将落在扇形内的点数m与所投点的总数n的比m/n作为k的近似值。P落在扇形内的充要条件是x^2+y^2<=1。

算法：

      1、确定产生点n的个数和缓冲区m（m<=n）的值，声明互斥锁

2、某一线程进入临界区，上锁

3、该线程一次性生成m个数，其他线程若想进入则挂起等待

4、该线程退出临界区，解锁 ，开始对刚才生成的随机点进行计算

5、重复2-4步，直至每个线程均完成对所要求点的操作

6、统计COUNTi的值

7、计算的值

程序代码：

//并行算法

 

#include <stdio.h>

#include <pthread.h>

#include <time.h>

#include <stdlib.h>

#include <sys/time.h>

#include <unistd.h>

 

 

long cs=0; //循环次数

long count=0; //主线程有效次数

long count_thread=0; //thread线程有效次数

struct timeval start, finish; //定义开始结束时间

double diffsec,diffusec;

long t;//每次生成数据数量

pthread_mutex_t mutex;

long double *data_thread,*data_main;

 

//void *thread(void *)

void thread(void)

{

   

   int i=0,j=0;

   double x,y;

   //long double *data_thread; 

   long double data_thread[t]; 

   for(i=0;i<cs/t;i++)

   {

   // data_thread=new (long double)[t];

       pthread_mutex_lock (&mutex);//lock

       for(j=0;j<t;j++)

           data_thread[j]=(long double)rand()/(long double)RAND_MAX;

       pthread_mutex_unlock(&mutex);//unlock

       for(j=0;j<t;j+=2)

       {

           x=data_thread[j];

           y=data_thread[j+1];

           if((x*x+y*y)<=1)

               count_thread++;

       }

       //delete []data_thread;

   }

}

 

//主线程计算量为总数的一半

int main(void)

{

   printf("Please input the number:");

   scanf("%d",&cs);

       cs=cs*1000000;

   printf("Please input the number for generating the data once:");

   scanf("%d",&t);

   t=t*100000;

   pthread_t id;

   int ret;

   srand( (unsigned)time( NULL ) );

   pthread_mutex_init(&mutex,NULL);//声明互斥锁

   ret=pthread_create(&id,NULL,(void *)thread,NULL); //创建thread线程

   gettimeofday(&start,NULL); //记录开始时间

   int i=0,j=0;

   double x,y;

   long double data_main[t];  

   for(i=0;i<cs/t;i++)

   {

   // data_main=new (long double)[t];

       pthread_mutex_lock(&mutex);//lock

       for(j=0;j<t;j++)

           data_main[j]=(long double)rand()/(long double)RAND_MAX;

       pthread_mutex_unlock(&mutex);//unlock

       for(j=0;j<t;j+=2)

       {

           x=data_main[j];

           y=data_main[j+1];

           if((x*x+y*y)<=1)

               count++;

       }

       //delete []data_main;

   }

   pthread_join(id,NULL); //两线程同步

   count+=count_thread;

   gettimeofday(&finish,NULL); //记录结束时间

   //计算所耗时间

   diffsec=(double)finish.tv_sec-start.tv_sec;

   diffusec=(double)(finish.tv_usec-start.tv_usec)/(double)1000000;

   diffsec+=diffusec;

   //输出结果

   printf("Cost time=%f second /n",diffsec);

   printf("roop times=%d /n",cs);

   printf("effective times=%d /n",count);

   printf("pi= %f /n",4*(double)count/(double)cs);

   return(0);

}