如何算精確圓週率 原網址http://forum.martinx.idv.tw/thread.php?threadid=575&boardid=49&styleid=1 __________________________________________________________    RAPPER 於 2002-02-02 07:33 在 舊源碼論壇上提出     救命阿~~~此提已經困擾我快半年了(快瘋了) 總是沒有頭緒 懇請幫幫忙  之前有寫一個類似的程式是1乘到1000(1000!) 使用10000進制來解  但元週率太複雜了 要分別算加減乘除 之前我用的級數 --> (1 - 1/3 + 1/5 - 1/7 +......0) * 4  就算跑到compiler爆炸也解不出準確值  所以重點是我不會寫此方程式----> pi / 4 = 4 * arcTan(1/5) - arcTan(1/239)  arcTan(x) = x - x^3 / 3 + x^5 / 5 - x^7 / 7 ......... --->要到第幾項才可停止???  整個程式大概要怎麼寫才較容易 ??????????????????  ___________________________________________________________________    macro 於 2002-02-02 07:33 在 舊源碼論壇上回應  ____________________________________________________________________ 應該用Taylor展開式就可以了吧  得出來的值與pi之間的誤差 也可以利用公式求得（可以用來檢驗精確度）  詳細的內容建議您看一下微積分課本（最好是原文的）  應該都有詳細的討論  瞭解Taylor展開式後  在程式裡加入一個累加各項的loop應該就可以了  當然結果的精確度與loop執行的次數成正比  這樣在compile的時候應該不會「爆掉」  倒是執行時若輸入的n值（精確度）過大  可能要算很久  比較有問題的應該是當顯示位數不足時  （int、long int、unsigned long int）  可能比較麻煩 要用到一點小技巧  不過應該也不難解決  如果要我寫這個程式的話  我可能會用 int array 來表示一個很長的整數吧     我其實沒寫過這樣的程式  是想試試 但是沒時間  在幾個討論區都有看到你的留言  所以小弟就現醜了  希望對你有點幫助  ______________________________________________________________________ parkus 於 2002-02-25 10:51 在 舊源碼論壇上回應     arcTan(x) = x - x^3 / 3 + x^5 / 5 - x^7 / 7 .........  =(x - x^3 / 3)+(x^5 / 5 - x^7 / 7)+.......  是一個交錯級數，每兩項決定誤差的大小。因此你要先設定你要誤差的大小(即精確度)，譬如因現在x=1/5及1/239，誤差的決定是由x=1/5 的項來決定。而當x=1/5，(x - x^3 / 3)=0.1973, (x^5 / 5 - x^7 / 7)=6.2171e-005, (x^9/9-x^11/11)=5.5027e-008,... 比較上面各項的大小，很明顯的， 你所採用的級數的最後兩項(被配對在一起的)的數量級(Order)就大約是誤差的大小。譬如採用三項的公式，arcTan(x) = (x - x^3 / 3) +(x^5 / 5 - x^7 / 7)+(x^9/9-x^11/11)，並不難看出所有忽略掉的後面項數的加總並不會比(x^9/9-x^11/11) 這項大，而這項的值大約是10^-8。所以小數點第8位就是你採用上面公式計算所得的arcTan(x)準確位數。 當然，從pi的公式，可得知大致這也是pi的精確位數。簡言之，若只要求算出pi的精確度在小數點8位，採用三項就可以了。若要十位或更多位精確度，四項以上就是必要的。     ps:  剛剛算了一下，使用上述配對的目前公式，算出的pi和正確的pi的數值如下：  pi=3.1415926526 (三項公式)  pi=3.1415926536 (真正值)  已經精確到第九位了。  若使用四項，  pi=3.1415926535886 (四項公式)  pi=3.1415926535898 (真正值)  則經精確到第十二位了 _____________________________________________________________________ RAPPER 於 2002-02-02 07:33 在 舊源碼論壇上回應     ****************************************** *以下是加州理工學院的學生寫的            * *可供大家參考喔                          * *算到小數 10000 位只要 10 秒!!!!!!!!!!!  * *以下是加州理工學院的學生寫的            * *可供大家參考喔                          * ****************************************** CODE:  /****************************************/  /* Compute pi to arbitrary precision    */  /* Author Roy Williams February 1994    */  /*   [EMAIL]roy@ccsf.caltech.edu[/EMAIL]               */  /* Uses Machin's formula...             */  /* pi/4 = 4 arctan(1/5) - arctan(1/239) */  /****************************************/  /* compile with cc -O -o pi1 pi1.c        */  /* run as "pi 1000" for 1000 digits     */  /****************************************/  /* The last few digits may be wrong.......... */  #include  #include #define BASE 10000 int nblock, *tot, *t1, *t2, *t3; void arctan(int *, int *, int *, int, int); void copy(int *, int *); bool zero(int *); void add(int *, int *); void sub(int *, int *); void mult(int *, int); void div(int *, int); void set(int *, int); void print(int *); void main(void) { int ndigit; printf("How many digit you want ?"); scanf("%d",&ndigit); if(ndigit < 20) ndigit = 20; nblock = ndigit/4; tot = (int *)malloc(nblock*sizeof(int)); t1 = (int *)malloc(nblock*sizeof(int)); t2 = (int *)malloc(nblock*sizeof(int)); t3 = (int *)malloc(nblock*sizeof(int)); if(!tot || !t1 || !t2 || !t3){ fprintf(stderr, "Not enough memory\n"); exit(1); } arctan(tot, t1, t2, 5, 1); mult(tot, 4); arctan(t3, t1, t2, 239, 2); sub(tot, t3); mult(tot, 4); print(tot); } void arctan(int *result, int *w1, int *w2, int denom, int onestep) { int denom2 = denom * denom; int k = 1; set(result, 1); div(result, denom); copy(w1, result); do{ if(onestep) div(w1, denom2); else { div(w1, denom); div(w1, denom); } copy(w2, w1); div(w2, 2*k 1); if(k%2) sub(result, w2); else add(result, w2); k ; } while(!zero(w2)); } void copy(int *result, int *from) { for(int i=0; i=0; i--){ result[i] = increm[i]; if(result[i] >= BASE){ result[i] -= BASE; result[i-1] ; } } } void sub(int *result, int *decrem) { int i; for(i=nblock-1; i>=0; i--){ result[i] -= decrem[i]; if(result[i] < 0){ result[i] = BASE; result[i-1]--; } } } void mult(int *result, int factor) { int i, carry = 0; for(i=nblock-1; i>=0; i--){ result[i] *= factor; result[i] = carry; carry = result[i]/BASE; result[i] %= BASE; } } void div(int *result, int denom) { int i, carry = 0; for(i=0; i

------
**********************************************************
哈哈&兵燹
最會的2大絕招 這個不會與那個也不會 哈哈哈 粉好

Delphi K.Top的K.Top分兩個字解釋Top代表尖端的意思,希望本討論區能提供Delphi的尖端新知
K.表Knowlege 知識,就是本站的標語:Open our mind