@[toc]

1 实验题目

计算π的一个有趣方法是使用一个称为Monte Carlo的技术，这种技术涉及随机。该技术工作如下：假设有一个圆，它内嵌一个正方形，如下图所示。

首先，通过(x,y)坐标生成一系列的随机点。这些点应在正方形内。在这些随机产生的点中，有的会落在圆内。接着，根据下面公式，估算π。
π= 4✖️（圆内的点数）/（总点数）

（1）请编写一个多线程版的算法，它创建一组子线程，每个子线程产生一组随机点，并确定点是否落在圆内。每个子线程更新所有落在园内的点的全局计数。当所有子线程计算完毕后，使用父线程计算π。（通过采用互斥锁，保护对共享全局变量的更新）？
（2）针对题目（1）请考虑一下，该怎样来提高计算π的性能（速度）？请实现相关的思想或算法的代码，并给出分析说明和实验数据证明！

2蒙特卡洛方法介绍

百度百科：百度百科
蒙特卡洛法，也称为统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。

3分析思路

①蒙特卡洛的本质是统计方法，所以必然要获取大量的数据才能得到较为精确的圆周率，
②在计算π的公式中，S=πrr，所以π=S/（rr），当r=1时，π此时等于S
③可以假设圆的半径为1，并且此时圆的外接正方形为22=4，设圆心位于坐标轴（0，0）④通过产生大量位[-1,1]的随机数坐标（x,y）,即-1<=x<=1,-1<=y<=1,若x^2+y^2<=1,则某个随机点落在圆内，否则落在圆外
⑤最终我们所求的圆周率π即为随机数的个数和/落在园内的随机数和。
⑥通过使用互斥锁来保证结果的正确性

4代码实现

4.1 安装插件

POSIX下开发多线程主要依赖的就是Pthread。使用它需要包含头文件#include<pthread.h>。因为这个库在Pthread之中，在编译的时候需要加上参数：-lpthread.

4.1.1

在终端输入命令

sudo apt-get install glibc-doc
然后输入密码

sudo apt-get install manpages-posix manpages-posix-dev
然后输入密码

自行百度linux下载vim（如果有就不用下载了）

4.1.2

新建一个文件夹tt
在文件夹中右击打开终端
输入

vi test2.c
按 i 进入编辑模式

4.2代码编写

#include<stdio.h>
#include<unistd.h>
#include<pthread.h>
#include<stdlib.h>
#include <time.h>

#define random_1(a,b) ((rand()%(b-a))+a)  //随机值将含a不含b
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
int value = 0;
long double valid=0,all=999999;//初始化有效点为0,总点数为999999
long double pi=0;
void*  fun(void* arg) {
    pthread_mutex_lock(&mutex);   
    long double x,y;//定义坐标x,y
    int i=1;
        srand((double)time(NULL)); //设置随机数种子，防止每次产生的随机数相同

    for(i=1;i<=999999;i++)
    {
    //获取x,y的随机值,0<=x<100000,0<=y<100000
    x=random_1(0,100000);
    y=random_1(0,100000);
    if(x*x+y*y<=10000000000)
    valid++;
    }
    pthread_mutex_unlock(&mutex);
    return (void*)0;    
}
//主函数
int main()
{
   clock_t start, end;
    start = clock();  //读取程序开始运行的时间
    pthread_t theradId1, threadId2;//定义两个子线程
    int error;//定义线程返回值error,若是成功建立线程返回0,否则返回错误的编号
    pthread_mutex_init(&mutex, NULL);//互斥锁初始化
    error = pthread_create(&theradId1, NULL, fun, NULL);
    sleep(2);//用sleep()控制子线程的不同并发过程
   
    if(error) {
        printf("Thread1 creat failed!\n");
        exit(1);
    }
   error = pthread_create(&threadId2, NULL, fun, NULL);
    if(error) {
        printf("Thread2 creat failed!\n");
        exit(1);
    }
   end = clock();           //记录结束时间
    pi=valid/999999*4;
    printf("valid的值为： %Lf\n",valid);
    printf("总的点数的值为999999\n");
    printf("π的值为： %Lf\n",pi);
    double seconds  =(double)(end - start)/CLOCKS_PER_SEC*1000;
     printf("Use time is: %.3fms\n", seconds);
    //进程锁销毁
    pthread_mutex_destroy(&mutex);
    return 0;
}

4.3执行代码

按 ESC 键 然后输入":“(冒号），再输入"wq”,回车

输入

gcc -o test test2.c -lpthread

再输入

./test

显示结果
当选取的点数较少是，例如200个，此时结果并不准确

当选取的总点数较多时，例如9999999个，结果相对更接近π

//可修改部分参数使结果更精确（多次测量求平均值）

5遇到的问题及解决方法

5.1随机数的获取问题

通过rand（）函数可以获取整数，如果半径控制在0~1的话，可能没办法实现随机小数获取，可以假设半径为100000（double型），当x,y的平方和小于10000000000时就可以说明在园内，否则在圆外

5.2随机数种子问题

如果不设置随机数种子，则每次获得的随机数在1s内相同，无法生成真正的随机数

6实验改进：OpenMP方法

6.1 OpenMP方法介绍

Open Multi-Processing的缩写，是一个应用程序接口（API），可用于显式指导多线程、共享内存的并行性。
在项目程序已经完成好的情况下不需要大幅度的修改源代码，只需要加上专用的pragma来指明自己的意图，由此编译器可以自动将程序进行并行化，并在必要之处加入同步互斥以及通信。当选择忽略这些pragma，或者编译器不支持OpenMp时，程序又可退化为通常的程序(一般为串行)，代码仍然可以正常运作，只是不能利用多线程来加速程序执行。OpenMP提供的这种对于并行描述的高层抽象降低了并行编程的难度和复杂度，这样程序员可以把更多的精力投入到并行算法本身，而非其具体实现细节。对基于数据分集的多线程程序设计，OpenMP是一个很好的选择。
OpenMP支持的语言包括C/C++、Fortran；而支持OpenMP的编译器VS、gcc、clang等都行。可移植性也很好：Unix/Linux和Windows

6.2 编程实现

打开终端输入

vi test2_2.c
按 i 进入编辑模式

输入代码

#include<stdio.h>
#include<unistd.h>
#include<pthread.h>
#include<stdlib.h>
#include <time.h>
#include<omp.h>
#define random_1(a,b) ((rand()%(b-a))+a)  //随机值将含a不含b
pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
int value = 0;
long double valid=0,all=999999;//初始化有效点为0,总点数为999999
long double pi=0;
//void*  fun(void* arg) {
void fun(){

   // pthread_mutex_lock(&mutex);   
    long double x,y;//定义坐标x,y

    int i=1;
              
        srand((double)time(NULL)); //设置随机数种子，防止每次产生的随机数相同

    for(i=1;i<=999999;i++)
    {


    //获取x,y的随机值,0<=x<100000,0<=y<100000
    x=random_1(0,100000);
    y=random_1(0,100000);
    if(x*x+y*y<=10000000000)
    #pragma omp crirical         //OpenMP并行编排支持
    {
    valid++;
    }
    }
   // pthread_mutex_unlock(&mutex);
       //  return (void*)0;    
}
//主函数
int main()
{
   clock_t start, end;
    start = clock();  //读取程序开始运行的时间

   fun();
   end = clock();           //记录结束时间
    pi=valid/999999*4;
    printf("valid的值为： %Lf\n",valid);
    printf("总的点数的值为999999\n");
    printf("π的值为： %Lf\n",pi);
    double seconds  =(double)(end - start)/CLOCKS_PER_SEC*1000;
     printf("Use time is: %.3fms\n", seconds);
    //进程锁销毁
    pthread_mutex_destroy(&mutex);
    return 0;
}

按Esc 再输入“：”，输入“wq”，回车
在终端输入

gcc -o test test2_2.c -fopenmp

再输入

./test

得到结果

和之前的普通互斥锁使用，在总点数相同的前提下（都为9999999个）可以明显看出使用OpenMP方法运行时间短很多

8参考内容

1,https://gist.github.com/lubobill1990/5388813
2,ubuntu如何安装pthread？
3,获取随机数
4,Linux多线程编程之pthread
5,如何输出long double？
6,linux shell 多线程执行程序
7,OpenMP入门教程1
8,linux下统计程序/函数运行时间

7源代码下载

链接密码: qsp3
仅供参考，请不要直接复制，还有很多不足之处，希望能共同进步，欢迎批评指正！