操作系统原理：C语言多线程加锁验证蒙特·卡罗（Monte Carlo）方法求π值

本文主要是介绍操作系统原理：C语言多线程加锁验证蒙特·卡罗（Monte Carlo）方法求π值，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

蒙特·卡罗方法

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是一类随机方法的统称。这类方法的特点是，可以在随机采样上计算得到近似结果，随着采样的增多，得到的结果是正确结果的概率逐渐加大。

在本实验中通过在正方形区域中生成随机点，记录随机点在圆形区域中的个数计算 $\pi$ 值。
$\pi= 4 \times (number \; of \; points \; in \; circle) / (total \; number \; of \; points)$
蒙特卡洛算法

整体思路

在主函数中开辟两个线程，在线程函数中分别将其绑定在两个CPU核上进行运算，通过一个循环源源不断地在正方形区域中生成随机点。

为了完成结果统计，定义了两个全局变量分别表示总点数和在圆形区域内的点数。

另外，为了避免两个线程同时对同一变量进行修改，需要对两个全局变量变量在使用的时候进行加锁。为了保证运行速度，在加锁的过程中要使临界区尽可能小。

【完整代码见文章最后】

生成随机点

在主函数中先使用srand(1)初始化随机数种子，在进程函数中通过对rand()函数的放缩和平移运算得到 -1~1 范围内的小数。

//生成随机点 [-1,1]
x = 2.0 * rand() / (double)RAND_MAX - 1;
y = 2.0 * rand() / (double)RAND_MAX - 1;

判断点在圆内

单独使用一个函数进行判断点是否在圆内，通过对传入的一组点坐标进行平方和运算得出判断结果。因为C语言中是没有 bool 类型的，所以这里“在圆内”则返回1，“不在圆内”则返回0 。

//判断点是否在圆内
int inCircle(double x, double y) {int flag = 1;if ((x * x + y * y) <= 1.0)flag = 1;elseflag = 0;return flag;
}

变量加锁

首先声明两个全局变量和分别与其对应的两个锁，并在主函数中将其初始化。

//全局变量
int sum_dots = 1000000;	//点的总数
int sum_in_circle = 0;	//在圆内的点的总数//声明锁
pthread_mutex_t lock_sum;	//对sum_dots的锁
pthread_mutex_t lock_in;	//对sum_in_circle的锁

//初始化锁
pthread_mutex_init(&lock_sum, NULL);
pthread_mutex_init(&lock_in, NULL);

在线程中有一个生成随机点的循环，每次循环开始的时候先将点总数sum_dots上锁pthread_mutex_lock(&lock_sum)，判断是否还有剩余点，如果有则将点数减一后立刻解锁pthread_mutex_unlock(&lock_sum);，如果没有剩余点也立刻解锁同时 break 出循环。

之后生成随机点，并判断点是否在圆中，若在圆中则对sum_in_circle上锁，进行加一运算后立刻解锁，最大程度上减小临界区域。

while (1){pthread_mutex_lock(&lock_sum);	//判断前，对点总数上锁if (sum_dots > 0) {sum_dots--;pthread_mutex_unlock(&lock_sum);//生成随机点 [-1,1]x = 2.0 * rand() / (double)RAND_MAX - 1;y = 2.0 * rand() / (double)RAND_MAX - 1;//如果生成点在圆中，圆内点总数+1if (inCircle(x, y)) {pthread_mutex_lock(&lock_in);sum_in_circle++;pthread_mutex_unlock(&lock_in);}}else {pthread_mutex_unlock(&lock_sum);break;}		}

结果截图

本实验中进行了 4 次小实验，分别为 100万个测试点时的单线程和双线程结果和 1000万个测试点时的单线程和双线程结果。

100万个点双线程：

100万个点单线程：

100万个点单线程

1000万个点双线程：

1000万个点双线程

1000万个点单线程：

1000万个点单线程

结果分析

当测试点个数相同时，从结果可以看出双线程的运行时间明显长于单线程。这是因为加锁缘故，其它线程在临界区内会有所暂停，导致了整体运行时间长于单线程。

至于同样测试点的条件下，双线程的精度高于单线程，可能是实验偶然性，但是我之后又做了几组实验同样是这样的结果。我猜测可能是因为C语言中的rand()函数是伪随机，加上我在实验中用的初始化随机种子是定值，因此单线程的结果可以复现，随机数据的混乱程度不高，但是多线程在调度rand随机序列的顺序上又多了一层随机性，可能提高了rand()函数的随机性。当然这个猜测并没有理论依据，希望路过的大佬给予解答。

当都是单线程或都是双线程时，测试点越多预测精度越高，误差越小，这个也可以从上面的结果中看出来，这即为概率论中的大数定律，也是蒙特·卡罗（Monte Carlo）方法的精髓所在。

代码

双线程用蒙特·卡罗（Monte Carlo）方法求 $\pi$ 值的完整代码如下，单线程方法仅在主函数中将其中一个线程注释掉即可。

#include <stdio.h>
#include <stdlib.h>
#ifndef __USE_GNU
#define __USE_GNU
#endif // !__USE_GNU
#include <unistd.h>
#include <sched.h>
#include <pthread.h>
#include <semaphore.h>double PI = 3.1415926535898;	//标准PI值//全局变量
int sum_dots = 1000000;	//点的总数
int sum_in_circle = 0;	//在圆内的点的总数//声明锁
pthread_mutex_t lock_sum;
pthread_mutex_t lock_in;//判断点是否在圆内
int inCircle(double x, double y) {int flag = 1;if ((x * x + y * y) <= 1.0)flag = 1;elseflag = 0;return flag;
}void* runner1() {	//将线程绑定到0号核上cpu_set_t cpuSet;CPU_ZERO(&cpuSet);CPU_SET(0, &cpuSet);sched_setaffinity(0, sizeof(cpuSet), &cpuSet);double x, y;	//随机点坐标while (1){pthread_mutex_lock(&lock_sum);	//判断前，对点总数上锁if (sum_dots > 0) {sum_dots--;pthread_mutex_unlock(&lock_sum);//生成随机点 [-1,1]x = 2.0 * rand() / (double)RAND_MAX - 1;y = 2.0 * rand() / (double)RAND_MAX - 1;//如果生成点在圆中，圆内点总数+1if (inCircle(x, y)) {pthread_mutex_lock(&lock_in);sum_in_circle++;pthread_mutex_unlock(&lock_in);}}else {pthread_mutex_unlock(&lock_sum);break;}		}pthread_exit(NULL);	//退出线程
}void* runner2() {//将线程绑定到1号核上cpu_set_t cpuSet;CPU_ZERO(&cpuSet);CPU_SET(1, &cpuSet);sched_setaffinity(0, sizeof(cpuSet), &cpuSet);double x, y;	//随机点坐标while (1) {pthread_mutex_lock(&lock_sum);	//判断前，对点总数上锁if (sum_dots > 0) {sum_dots--;pthread_mutex_unlock(&lock_sum);//生成随机点 [-1,1]x = 2.0 * rand() / (double)RAND_MAX - 1;y = 2.0 * rand() / (double)RAND_MAX - 1;//如果生成点在圆中，圆内点总数+1if (inCircle(x, y)) {pthread_mutex_lock(&lock_in);sum_in_circle++;pthread_mutex_unlock(&lock_in);}}else {pthread_mutex_unlock(&lock_sum);break;}}pthread_exit(NULL);	//退出线程
}int main(){int sum_dots_p = sum_dots;	//复制总点数，作最后计算用pthread_t tid1, tid2;		//线程IDpthread_attr_t attr;		//线程属性pthread_attr_init(&attr);	//设置默认线程属性//初始化随机数发生器 srand(1);//初始化锁pthread_mutex_init(&lock_sum, NULL);pthread_mutex_init(&lock_in, NULL);//执行两个线程分别进行随机生成点pthread_create(&tid1, &attr, runner1, NULL);pthread_create(&tid2, &attr, runner2, NULL);//等待两个线程pthread_join(tid1, NULL);pthread_join(tid2, NULL);//计算结果double estimate_PI = (double)(4.0 * sum_in_circle / sum_dots_p);printf("PI: %lf\n", estimate_PI);printf("Error Value: %lf\n", estimate_PI - PI);return 0;
}