蒙特卡洛法，多线程求pi 当启用多线程时，pi会偏小，可能因为线程切换对随机数生成有影响

README.md

计算PI的方法 蒙特卡洛法

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

蒙特卡洛法的工作思想

1. 用蒙地卡羅方法模拟某一过程时，需要产生各种概率分布的随机变量
2. 用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解

蒙特卡洛方法计算PI

蒙特卡洛方法的本质是统计模拟，在计算PI中，根据圆面积计算公式S=PI*r^2，得到PI的表达式为PI=S/r^2。在模拟过程如下

1. 随机生成[0,1]范围的一对随机数
2. 计算该点到(0,0)坐标的欧氏距离
3. 如果距离小于1，则表示落入圆内
4. 重复n次，得到点落入圆内的次数为m
5. PI的计算公式为 PI=4*m/n

程序代码主要部分

double e_distance(double a,double b){
    return sqrt((a-0.5)*(a-0.5)+(b-0.5)*(b-0.5));
}

void * calculate_pi(void * threadarg){
    srand((unsigned)time(0));
    int index;
    int calc_sum=0;
    int hit=0;
    double a,b;
    for (index = 0; index <loopNum; index++) {
        ++calc_sum;
        a=rand()/(double)max_int;
        b=rand()/(double)max_int;
        if(e_distance(a,b)<0.500000){
            ++hit;
        }  
    }  
    double res= (double)hit/calc_sum*4;
    resArr[cur_index++]=res;
}

执行结果

在程序中，我设置循环次数为100000000 得出的结果为:3.1415558

布丰投针法

18世纪，布丰提出以下问题：设我们有一个以平行且等距木纹铺成的地板（如右图），现在随意抛一支长度比木纹之间距离小的针，求针和其中一条木纹相交的概率。这就是布丰投针问题。 实际上，布丰投针法也应该算作一种蒙特卡洛方法，都是通过统计模拟的方法计算一个难以计算的值。

计算PI的步骤

1. 确定在二维平面上等距直线的距离为d
2. 随机生成一根针，针的长度为a，角度为b，位置为c
3. 计算这根针是否落在直线上
4. 重复n次，得出落在直线上的次数为m，针的总长度为l
5. 计算可得，PI可以表示为2l/(d*(m/n))

程序主要部分

void needle_range(double length,double angle,double position,double &top,double&bottom){
    top=abs(sin(angle)*length)+position;
    bottom=position;
}
/**
 * 生成两个随机数a,b
 * 计算和0.5,0.5的距离
 * 每次将calc_sum++
 * 如果小于0.5，则将hit++
 * pi*r*r/(4*r*r)*4
 */
void * calculate_pi(void * threadarg){
    srand((unsigned)time(0));
    int index;
    int calc_sum=0;
    int hit=0;
    double a,b,c;
    //a表示针的长度，介于0-1之间
    //b表示角度，是一个较大的double随机数
    //c表示位置，是一个较大的double随机数
    double top,bottom;
    double sum_length=0;
    for (index = 0; index <loopNum; index++) {
        ++calc_sum;

        a=rand()/(double)max_int;
        b=rand()/(double)10000;
        c=rand()/(double)10000;
        needle_range(a,b,c,top,bottom);
        if(ceil(bottom)==floor(top)){
            ++hit;
        }   
        sum_length+=a;

    }   
    double res= (double)2*sum_length/calc_sum/(hit/(double)calc_sum);
    resArr[cur_index++]=res;
}

执行结果

在程序中，循环次数为100000000 得出的结果为:3.14141051

随机三角形锐角概率法

随机生成能组成三角形的三条边，计算其组成锐角三角形的概率。

步骤

1. 在二维空间随机生成三个点p1,p2,p3
2. 计算三个点之间的距离,从大到小分别为a,b,c
3. 如果aa<bb+c*c，随机生成的三个点即可组成锐角三角形
4. 循环很多次，求得锐角三角形的概率为p
5. PI即可表示为PI=4-2*P

主要代码

/**
 * 根据三条边长确定能否组成锐角三角形
 */
bool compose_triangle(double a,double b,double c,bool &isAcute){
    double tmp;
    if(a<b){
        swap(a,b);
    }
    if(a<c){
        swap(a,c);
    }
    isAcute=(a*a)<(b*b+c*c);
    return (a<b+c);
}

void * calculate_pi(void * threadarg){
    srand((unsigned)time(0));
    int index;
    int calc_sum=0;
    int hit=0;
    double a,b,c;
    bool isAcute;
    for (index = 0; index <loopNum; ) {
        a=rand()/(double)max_int;
        b=rand()/(double)max_int;
        c=rand()/(double)max_int;

        if(compose_triangle(a,b,c,isAcute)){
            ++calc_sum;
            ++index;
            if(isAcute){
                ++hit;
            }
        }
    }
    double res=4-2*hit/(double)calc_sum;
    resArr[cur_index++]=res;
}

执行结果

在程序中，主要代码的循环次数为:10000000 执行结果为:3.1415968

buffon_pi.cpp

#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
const int max_int=std::numeric_limits<int>::max();
const int loopNum=100000000;
const int threadNum=1;
double resArr[threadNum];
int cur_index=0;
void *do_nothing(void *null) {
    int i;
    i=0;
    pthread_exit(NULL);
}                      
double e_distance(double a,double b){
    return sqrt((a-0.5)*(a-0.5)+(b-0.5)*(b-0.5));
}
void needle_range(double length,double angle,double position,double &top,double&bottom){
    top=abs(sin(angle)*length)+position;
    bottom=position;
}
/**
 * 生成两个随机数a,b
 * 计算和0.5,0.5的距离
 * 每次将calc_sum++
 * 如果小于0.5，则将hit++
 * pi*r*r/(4*r*r)*4
 */
void * calculate_pi(void * threadarg){
    srand((unsigned)time(0));
    int index;
    int calc_sum=0;
    int hit=0;
    double a,b,c;
    //a表示针的长度，介于0-1之间
    //b表示角度，是一个double随机数
    //c表示位置，是一个double随机数
    double top,bottom;
    double sum_length=0;
    for (index = 0; index <loopNum; index++) {
        ++calc_sum;

        a=rand()/(double)max_int;
        b=rand()/(double)10000;
        c=rand()/(double)10000;
        needle_range(a,b,c,top,bottom);
        if(ceil(bottom)==floor(top)){
            ++hit;
        }
        sum_length+=a;

    }
    double res= (double)2*sum_length/calc_sum/(hit/(double)calc_sum);
    resArr[cur_index++]=res;
}
void printResult(){
    double res=0;
    int i;
    for (i = 0; i < threadNum; i++) {
        res+=resArr[i];
    }
    cout<<"After "<<loopNum<<" loops, we calculate the pi is: "<<setprecision(10)<<res/threadNum<<endl;

}
int main(int argc, char *argv[]) {

    //struct thread_data thread_data_array[4];
    pthread_t threads[threadNum];
    int b=1;
    int * a=&b;
    int rc;
    int i;
    pthread_attr_t attr;
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

    for (i = 0; i < threadNum; i++) {
        rc=pthread_create(&threads[i],NULL,calculate_pi,a);
        if (rc) {
            printf("ERROR; return code from pthread_create() is %d\n", rc);
            exit(-1);
        }
    }
    pthread_attr_destroy(&attr);
    void* status=malloc(sizeof(long));
    for (i = 0; i < threadNum; i++) {
        rc=pthread_join(threads[i],&status);
        if (rc) {
            printf("ERROR; return code from pthread_join() is %d\n", rc);
            exit(-1);
        }
        printf("Main: completed join with thread %d having a status   of %ld\n",i,(long)status);
    }
    printf("Main: program completed. Exiting.\n");
    printResult();
    pthread_exit(NULL);
}

montecarlo_pi.cpp

#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
const int max_int=std::numeric_limits<int>::max();
const int loopNum=100000000;
const int threadNum=1;
double resArr[threadNum];
int cur_index=0;
void *do_nothing(void *null) {
    int i;
    i=0;
    pthread_exit(NULL);
}                      
double e_distance(double a,double b){
    return sqrt((a-0.5)*(a-0.5)+(b-0.5)*(b-0.5));
}
/**
 * 生成两个随机数a,b
 * 计算和0.5,0.5的距离
 * 每次将calc_sum++
 * 如果小于0.5，则将hit++
 * pi*r*r/(4*r*r)*4
 */
void * calculate_pi(void * threadarg){
    srand((unsigned)time(0));
    int index;
    int calc_sum=0;
    int hit=0;
    double a,b;
    for (index = 0; index <loopNum; index++) {
        ++calc_sum;
        a=rand()/(double)max_int;
        b=rand()/(double)max_int;
//        cout<<a<<"\t"<<b<<"\t"<<distance<<endl;
        if(e_distance(a,b)<0.500000){
            ++hit;
        }
    }
    double res= (double)hit/calc_sum*4;
    resArr[cur_index++]=res;
}
void printResult(){
    double res=0;
    int i;
    for (i = 0; i < threadNum; i++) {
        res+=resArr[i];
    }
    cout<<"After "<<loopNum<<" loops, we calculate the pi is: "<<setprecision(10)<<res/threadNum<<endl;
 
}
int main(int argc, char *argv[]) {
    //struct thread_data thread_data_array[4];
    pthread_t threads[threadNum];
    int b=1;
    int * a=&b;
    int rc;
    int i;
    pthread_attr_t attr;
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

    for (i = 0; i < threadNum; i++) {
        rc=pthread_create(&threads[i],NULL,calculate_pi,a);
        if (rc) {
            printf("ERROR; return code from pthread_create() is %d\n", rc);
            exit(-1);
        }
    }
    pthread_attr_destroy(&attr);
    void* status=malloc(sizeof(long));
    for (i = 0; i < threadNum; i++) {
        rc=pthread_join(threads[i],&status);
        if (rc) {
            printf("ERROR; return code from pthread_join() is %d\n", rc);
            exit(-1);
        }
        printf("Main: completed join with thread %d having a status   of %ld\n",i,(long)status);
    }
    printf("Main: program completed. Exiting.\n");
    printResult();
    pthread_exit(NULL);
}

triangle_pi.cpp

#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
const int max_int=std::numeric_limits<int>::max();
const int loopNum=10000000;
const int threadNum=1;
double resArr[threadNum];
int cur_index=0;
void *do_nothing(void *null) {
    int i;
    i=0;
    pthread_exit(NULL);
}                      
double e_distance(double a,double b){
    return sqrt((a-0.5)*(a-0.5)+(b-0.5)*(b-0.5));
}
void needle_range(double length,double angle,double position,double &top,double&bottom){
    top=abs(sin(angle)*length)+position;
    bottom=position;
}
inline void swap(double *a,double *b){
    double tmp;
    tmp=*a;
    *a=*b;
    *b=tmp;
}
/**
 * 根据三条边长确定能否组成锐角三角形
 */
bool compose_triangle(double a,double b,double c,bool &isAcute){
    double tmp;
    if(a<b){
        swap(a,b);
    }
    if(a<c){
        swap(a,c);
    }
    isAcute=(a*a)<(b*b+c*c);
    return (a<b+c);
}
/**
 * 生成两个随机数a,b
 * 计算和0.5,0.5的距离
 * 每次将calc_sum++
 * 如果小于0.5，则将hit++
 * pi*r*r/(4*r*r)*4
 */
void * calculate_pi(void * threadarg){
    srand((unsigned)time(0));
    int index;
    int calc_sum=0;
    int hit=0;
    double a,b,c;
    bool isAcute;
    for (index = 0; index <loopNum; ) {
        a=rand()/(double)max_int;
        b=rand()/(double)max_int;
        c=rand()/(double)max_int;

        if(compose_triangle(a,b,c,isAcute)){
            ++calc_sum;
            ++index;
    //        cout<<a<<"\t"<<b<<"\t"<<c<<endl;
            if(isAcute){
                ++hit;
            }
        }
    }
    double res=4-2*hit/(double)calc_sum;
    resArr[cur_index++]=res;
}
void printResult(){
    double res=0;
    int i;
    for (i = 0; i < threadNum; i++) {
        res+=resArr[i];
    }
    cout<<"After "<<loopNum<<" loops, we calculate the pi is: "<<setprecision(10)<<res/threadNum<<endl;

}
int main(int argc, char *argv[]) {

    //struct thread_data thread_data_array[4];
    pthread_t threads[threadNum];
    int b=1;
    int * a=&b;
    int rc;
    int i;
    pthread_attr_t attr;
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

    for (i = 0; i < threadNum; i++) {
        rc=pthread_create(&threads[i],NULL,calculate_pi,a);
        if (rc) {
            printf("ERROR; return code from pthread_create() is %d\n", rc);
            exit(-1);
        }
    }
    pthread_attr_destroy(&attr);
    void* status=malloc(sizeof(long));
    for (i = 0; i < threadNum; i++) {
        rc=pthread_join(threads[i],&status);
        if (rc) {
            printf("ERROR; return code from pthread_join() is %d\n", rc);
            exit(-1);
        }
        printf("Main: completed join with thread %d having a status   of %ld\n",i,(long)status);
    }
    printf("Main: program completed. Exiting.\n");
    printResult();
    pthread_exit(NULL);
}