藍影

前言

對於數學家和程式設計者，計算 PI 的準確度一直都是值得挑戰的問題。我在此提出三個正常求PI的方式，及另一個覺得非常 Kuso 的求 PI 方式。再次強調，求 PI 的方法真的非常多，同時 paper 也有不斷的更新，我提供的只是之前學者所研究出來的方法，實際上要算小數點後幾百、幾千、幾萬位數的 PI，請配合其它的演算法進行。 

最後將附上程式碼，再與 google 上所附的 PI 值做比較。

歐拉(Euler)等式求PI

 當然，這個公式用 C 語言比較麻煩的是還要調用 sqrt() 函數。此外類似這種公式的型式實在是多到不行，比如說：

至於其它的公式我就不再此贅述了..
但使用 Euler 等式求PI，其收斂性不佳。

Leibniz 求PI
  

Leibniz 的求法並不能算是正確，只能算是 "精簡"，其準確度仍待商議，
同時收斂性和 Euler 一樣，也是不佳。

Machin求PI

這個式子應該是最簡單求 PI 的方式，只有一個式子，用不到 Loop，同時精準度也高。

Kuso - Monte_Carlo 求 PI

蒙地卡羅法求 PI，是個非常賭博性的求法。
隨機抓一個點 (x, y)，判斷 (x,y) 是不是在圓內，並隨機測試 N 次，在圓內的次數為 n 次，則 PI = 4 * (n/N)。

程式碼

// ====================================
// FileName: PI.cpp
// Author  : Edison.Shih.
// Complier: VC 2008

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h> // use Monte_Carlo_PI

#define PI (double)(3.14159265)

// ====================================
// declare function

double Euler_PI(unsigned long n);
double Leibniz_PI(unsigned long n);
double Machin_PI();
double Monte_Carlo_PI(unsigned long n);

// ====================================
// Euler: PI*PI / 6 = 1/1*1 + 1/2*2 + 1/3*3 +.......

double Euler_PI(unsigned long n){
        unsigned long i=0;
        double ret = 0.0;
        for(i=1; i<=n; i++) ret += 1.0 / (double)(i*i);
        return sqrt(ret*6.0);
}

// ====================================
// Leibniz: PI = 4* (1/1 - 1/3 + 1/5 - 1/7 +....)
double Leibniz_PI(unsigned long n)
{
        unsigned long i=0;
        double ret = 0.0;
        for(i=0; i<n; i++){
                if(!(i%2)) ret = ret + 1.0/(double)(2*i+1);
                else ret = ret - 1.0/(double)(2*i+1);
        }
        return (4.0*ret);
}

// ====================================
// Machin: PI = 16*arctan(1/5) - 4*arctan(1/239)

double Machin_PI(){
        return (16*atan(1.0/5.0) - 4*atan(1.0 / 239.0));
}

// ====================================
// Monte_Carlo_PI
double Monte_Carlo_PI(unsigned long n){
        const unsigned MAX = 10000;
        unsigned long i=0;
        double x=0.0, y=0.0;
        unsigned long in_times = 0;

        srand(time(NULL));
        for(i=0; i<n ;i++){
                x = rand() % MAX;
                y = rand() % MAX;
                if(x*x + y*y < (double)(MAX*MAX)) in_times++;
        }
        return (4.0 * (double)(in_times)/(double)(n));
}

// ====================================
// main function
int main(int argc, char **argv)
{
        printf("PI = %lf\n", PI);
        printf("Euler_PI(100) = %lf\n", Euler_PI(100)); // 展開100項
        printf("Leibniz_PI(100) = %lf\n", Leibniz_PI(100)); // 展開100項
       
        printf("Machin_PI = %lf\n", Machin_PI());
        printf("Monte_Carlo_PI(100000) = %lf\n", Monte_Carlo_PI(100000)); // 一萬次隨機測試

}
 // ====================================

輸出結果

PI = 3.141593
Euler_PI(100) = 3.132077
Leibniz_PI(100) = 3.131593
Machin_PI = 3.141593
Monte_Carlo_PI(100000) = 3.279160

後序

結論很明顯看得出來，使用Machin 的方法較簡單，而且也較其它方法準確、快！

至於Monte_Carlo的方法...嗯...那就留給各位評論了...

最後聲明，我不是數學家，公式的推導過程我並不清楚，但這些公式在網路上的確不難找到，我只是將它試著用C語言實剛已作出來而已，如果有問題，也歡迎討論一起討論。