圓周率后的小數位數是無止境的，如何使用電腦來計算這無止境的小數是一些數學家與程式設計師所感興趣的，在這邊介紹一個公式配合 大數運算，可以計算指定位數的圓周率。

解法首先介紹J.Marchin的圓周率公式：

PI = [16/5 - 16 / (3*53) + 16 / (5*55) - 16 / (7*57) + ......] -

??????[4/239 - 4/(3*2393) + 4/(5*2395) - 4/(7*2397) + ......]

?

可以將這個公式整理為：

PI = [16/5 - 4/239] - [16/(53) - 4/(2393)]/3+ [16/(55) - 4/(2395)]/5 + ......

?

也就是說第n項，若為奇數則為正數，為偶數則為負數，而項數表示方式為：

[16/52*n-1?- 4/2392*n-1] / (2*n-1)

?

如果我們要計算圓周率至10的負L次方，由于[16/52*n-1?- 4/2392*n-1]中16/52*n-1比4/2392*n-1來的大，具有決定性，所以表示至少必須計算至第n項：

[16/52*n-1?] / (2*n-1) = 10-L

?

將上面的等式取log并經過化簡，我們可以求得：

n = L / (2log5) = L / 1.39794

?

所以若要求精確度至小數后L位數，則只要求至公式的第n項，其中n等于：

n = [L/1.39794] + 1

?

在上式中[]為高斯符號，也就是取至整數（不大于L/1.39794的整數）；為了計簡方便，可以在程式中使用下面這個公式來計簡第n項：

[Wn-1/52- Vn-1 / (2392)] / (2*n-1)

?

這個公式的演算法配合大數運算函式的演算法為： div(w, 25, w);

div(v, 239, v);

div(v, 239, v);

sub(w, v, q);

div(q, 2*k-1, q)

?

至于大數運算的演算法，請參考之前的文章，必須注意的是在輸出時，由于是輸出陣列中的整數值，如果陣列中整數位數不滿四位，則必須補上0，在C語言中只要 使用格式指定字%04d，使得不足位數部份自動補上0再輸出，至于Java的部份，使用 NumberFormat來作格式化。

#include <stdio.h> #define L 1000 #define N L/4+1 // L 為位數，N是array長度 void add(int*, int*, int*); void sub(int*, int*, int*); void div(int*, int, int*); int main(void) { int s[N+3] = {0}; int w[N+3] = {0}; int v[N+3] = {0}; int q[N+3] = {0}; int n = (int)(L/1.39793 + 1); int k; w[0] = 16*5; v[0] = 4*239; for(k = 1; k <= n; k++) { // 套用公式 div(w, 25, w); div(v, 239, v); div(v, 239, v); sub(w, v, q); div(q, 2*k-1, q); if(k%2) // 奇數項 add(s, q, s); else // 偶數項 sub(s, q, s); } printf("%d.", s[0]); for(k = 1; k < N; k++) printf("%04d", s[k]); printf("\n"); return 0; } void add(int *a, int *b, int *c) { int i, carry = 0; for(i = N+1; i >= 0; i--) { c[i] = a[i] + b[i] + carry; if(c[i] < 10000) carry = 0; else { // 進位 c[i] = c[i] - 10000; carry = 1; } } } void sub(int *a, int *b, int *c) { int i, borrow = 0; for(i = N+1; i >= 0; i--) { c[i] = a[i] - b[i] - borrow; if(c[i] >= 0) borrow = 0; else { // 借位 c[i] = c[i] + 10000; borrow = 1; } } } void div(int *a, int b, int *c) { // b 為除數 int i, tmp, remain = 0; for(i = 0; i <= N+1; i++) { tmp = a[i] + remain; c[i] = tmp / b; remain = (tmp % b) * 10000; } }

?

《新程序員》：云原生和全面數字化實踐50位技術專家共同創作，文字、視頻、音頻交互閱讀

總結

以上是生活随笔為你收集整理的长 PI的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。