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

解法首先介紹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項，若為奇數(shù)則為正數(shù)，為偶數(shù)則為負(fù)數(shù)，而項數(shù)表示方式為：

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

如果我們要計算圓周率至10的負(fù)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并經(jīng)過化簡，我們可以求得：

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

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

n = [L/1.39794] + 1

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

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

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

div(v, 239, v);

div(v, 239, v);

sub(w, v, q);

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

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

代碼部分

#include <stdio.h> #define L 1000 #define N L/4+1 // L 為位數(shù)，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) // 奇數(shù)項 add(s, q, s); else // 偶數(shù)項 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 為除數(shù) int i, tmp, remain = 0; for(i = 0; i <= N+1; i++) { tmp = a[i] + remain; c[i] = tmp / b; remain = (tmp % b) * 10000; } }

總結(jié)

以上是生活随笔為你收集整理的计算指定位数的圆周率的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。