linux環境下使用Monte Carlo計算π

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之中，在編譯的時候需要加上參數：-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 #include #include #include #include #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 #include #include #include #include #include #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

僅供參考，請不要直接復制，還有很多不足之處，希望能共同進步，歡迎批評指正！

Linux 任務調度

版權聲明：本文內容由網絡用戶投稿，版權歸原作者所有，本站不擁有其著作權，亦不承擔相應法律責任。如果您發現本站中有涉嫌抄襲或描述失實的內容，請聯系我們jiasou666@gmail.com 處理，核實后本網站將在24小時內刪除侵權內容。

版權聲明：本文內容由網絡用戶投稿，版權歸原作者所有，本站不擁有其著作權，亦不承擔相應法律責任。如果您發現本站中有涉嫌抄襲或描述失實的內容，請聯系我們jiasou666@gmail.com 處理，核實后本網站將在24小時內刪除侵權內容。