蒙特卡洛方法實現計算圓周率的方法比較簡單，其思想是假設我們向一個正方形的標靶上隨機投擲飛鏢，靶心在正中央，標靶的長和寬都是2 英尺。同時假設有一個圓與標靶內切。圓的半徑是1英尺，面積是π平方英尺。如果擊中點在標靶上是均勻分布的（我們總會擊中正方形），那么飛鏢擊中圓的數量近似滿足等式
飛鏢落在圓內的次數/飛鏢落在標靶內的總次數=π/4
因為環包含的面積與正方形面積的比值是π/4。
因為環所包含的面積與正方形面積的比值是π/4。
我們可以用這個公式和隨機數產生器來估計π的值。
偽代碼如下：
number_in_circle=0; for(toss=0;toss<number_of_tosses;toss++){ x=random double between -1 and 1; y=random double between -1 and 1; distance_squared=x*x+y*y; if(distance_squared<=1) number_in_cycle++; } pi_estimate=4*number_in_cycle/((double)number_in_tosses);
這種采用了隨機（隨機投擲）的方法稱為蒙特卡洛（Monte Carlo）方法。
編寫了一個采用蒙特卡洛方法的MPI，Pthread，OpenMP程序估計π的值。
使用MPI編寫時，進程0讀取總的投擲次數，並把它們廣播給各個進程。使用MPI_Reduce求出局部變量number_in_cycle的全局總和，並讓進程0打印它。
使用Pthread編寫時，有主線程讀入總的投擲數，然后輸出估算值。
使用OpenMP編寫時，在開啟任何線程前讀取總的投擲次數。使用reduction子句計算飛鏢集中環內的次數。在合並所有的線程后，打印結果。
這三個程序中，投擲次數可能都非常大，可能總的投擲次數和擊中環內的次數都得用long int型來表示。
下面是MPI程序代碼
1 #include<stdio.h> 2 #include<stdlib.h> 3 #include<math.h> 4 #include<time.h> 5 #include<mpi.h> 6 void read_num(long long int *num_point,int my_rank,MPI_Comm comm); 7 void compute_pi(long long int num_point,long long int* num_in_cycle,long long int* local_num_point,int comm_sz,long long int *total_num_in_cycle,MPI_Comm comm,int my_rank); 8 int main(int argc,char** argv){ 9 long long int num_in_cycle,num_point,total_num_in_cycle,local_num_point; 10 int my_rank,comm_sz; 11 MPI_Comm comm; 12 MPI_Init(NULL,NULL);//初始化 13 comm=MPI_COMM_WORLD; 14 MPI_Comm_size(comm,&comm_sz);//得到進程總數 15 MPI_Comm_rank(comm,&my_rank);//得到進程編號 16 read_num(&num_point,my_rank,comm);//讀取輸入數據 17 compute_pi(num_point,&num_in_cycle,&local_num_point,comm_sz,&total_num_in_cycle,comm,my_rank); 18 MPI_Finalize(); 19 return 0; 20 } 21 void read_num(long long int* num_point,int my_rank,MPI_Comm comm){ 22 if(my_rank==0){ 23 printf("please input num in sqaure \n"); 24 scanf("%lld",num_point); 25 } 26 /* 27 廣播函數 28 int MPI_Bcast( 29 void* data_p //in/out 30 int count //in 31 MPI_Datatype datatype //in 32 int source_proc //in 33 MPI_Comm comm //in 34 ) 35 */ 36 MPI_Bcast(num_point,1,MPI_LONG_LONG,0,comm); 37 38 } 39 void compute_pi(long long int num_point,long long int* num_in_cycle,long long int* local_num_point,int comm_sz,long long int *total_num_in_cycle,MPI_Comm comm,int my_rank){ 40 *num_in_cycle=0; 41 *local_num_point=num_point/comm_sz; 42 double x,y,distance_squared; 43 srand(time(NULL)); 44 for(long long int i=0;i< *local_num_point;i++){ 45 x=(double)rand()/(double)RAND_MAX; 46 x=x*2-1; 47 y=(double)rand()/(double)RAND_MAX; 48 y=y*2-1; 49 distance_squared=x*x+y*y; 50 if(distance_squared<=1) 51 *num_in_cycle=*num_in_cycle+1; 52 } 53 /* 54 全局函數 55 MPI_Reduce( 56 void* input_data_p //in 57 void* output_data_p //out 58 int count //in 59 MPI_Datatype datatype //in 60 MPI_Op oprtator //in 61 int dest_process //in 62 MPI_Comm comm //in 63 ) 64 */ 65 MPI_Reduce(num_in_cycle,total_num_in_cycle,1,MPI_LONG_LONG,MPI_SUM,0,comm); 66 if(my_rank==0){ 67 double pi=(double)*total_num_in_cycle/(double)num_point*4; 68 printf("the estimate value of pi is %lf\n",pi); 69 } 70 }
進行編譯 mpicc -std=c99 -o mpi_mete mpi_mete.c
執行 mpiexec -n 4 mpi_mete
輸入數據
下面是Pthread程序代碼
1 #include<stdlib.h> 2 #include<stdio.h> 3 #include<math.h> 4 #include<pthread.h> 5 int thread_count; 6 long long int num_in_circle,n; 7 pthread_mutex_t mutex; 8 void* compute_pi(void* rank); 9 int main(int argc,char* argv[]){ 10 long thread; 11 pthread_t* thread_handles; 12 thread_count=strtol(argv[1],NULL,10); 13 printf("please input the number of point\n"); 14 scanf("%lld",&n); 15 thread_handles=(pthread_t*)malloc(thread_count*sizeof(pthread_t)); 16 pthread_mutex_init(&mutex,NULL); 17 for(thread=0;thread<thread_count;thread++){ 18 //創建線程 19 /* 20 int pthread_create( 21 pthread_t* thread_p //out 22 const pthread_attr_t* attr_p 23 void* (*start_routine)(void*) //in 24 void* arg_p //in 25 ) 26 第一個參數是一個指針，指向對應的pthread_t對象。注意，pthread_t對象不是pthread_create函數分配的，必須在調用pthread_create函數前就為pthread_t 27 對象分配內存空間。第二個參數不用，所以只是函數調用時把NULL傳遞參數。第三個參數表示該線程將要運行的函數；最后一個參數也是一個指針，指向傳給函數start_routine的參數 28 */ 29 pthread_create(&thread_handles[thread],NULL,compute_pi,(void*)thread); 30 } 31 //停止線程 32 /* 33 int pthread_join( 34 pthread_t thread /in 35 void** ret_val_p /out 36 ） 37 第二個參數可以接收任意由pthread_t對象所關聯的那個線程產生的返回值。 38 */ 39 for(thread=0;thread<thread_count;thread++){ 40 pthread_join(thread_handles[thread],NULL); 41 } 42 pthread_mutex_destroy(&mutex); 43 double pi=4*(double)num_in_circle/(double)n; 44 printf("the esitimate value of pi is %lf\n",pi); 45 } 46 void* compute_pi(void* rank){ 47 long long int local_n; 48 local_n=n/thread_count; 49 double x,y,distance_squared; 50 for(long long int i=0;i<local_n;i++){ 51 x=(double)rand()/(double)RAND_MAX; 52 y=(double)rand()/(double)RAND_MAX; 53 distance_squared=x*x+y*y; 54 if(distance_squared<=1){ 55 pthread_mutex_lock(&mutex); 56 num_in_circle++; 57 pthread_mutex_unlock(&mutex); 58 } 59 } 60 return NULL; 61 }
進行編譯 gcc -std=c99 -o pth_mete pth_mete.c
運行 ./pth_mete 4
輸入數據結果如下
下面是OpenMP代碼
1 #include<stdlib.h> 2 #include<stdio.h> 3 #include<time.h> 4 #include<omp.h> 5 int main(int argc,char** argv){ 6 long long int num_in_cycle=0; 7 long long int num_point; 8 int thread_count; 9 thread_count=strtol(argv[1],NULL,10); 10 printf("please input the number of point\n"); 11 scanf("%lld",&num_point); 12 srand(time(NULL)); 13 double x,y,distance_point; 14 long long int i; 15 #pragma omp parallel for num_threads(thread_count) default(none) \ 16 reduction(+:num_in_cycle) shared(num_point) private(i,x,y,distance_point) 17 for( i=0;i<num_point;i++){ 18 x=(double)rand()/(double)RAND_MAX; 19 y=(double)rand()/(double)RAND_MAX; 20 distance_point=x*x+y*y; 21 if(distance_point<=1){ 22 num_in_cycle++; 23 } 24 } 25 double estimate_pi=(double)num_in_cycle/num_point*4; 26 printf("the estimate value of pi is %lf\n",estimate_pi); 27 return 0; 28 }
編譯代碼 gcc -std=c99 -fopenmp -o omp_mete omp_mete.c
執行 ./omp_mete 4
結果如下
