圓周率後的小數位數是無止境的,如何使用電腦來計算這無止境的小數是一些數學家與程式設計師所感興趣的,在這邊介紹一個公式配合 大數運算,可以計算指定位數的圓周率。
解法首先介紹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);
sub(w, v, q);
div(q, 2*k-1, q)
至于大數運算的演算法,請參考之前的文章,必須注意的是在輸出時,由于是輸出陣列中的整數值,如果陣列中整數位數不滿四位,則必須補上0,在C語言中隻要 使用格式指定字%04d,使得不足位數部份自動補上0再輸出,至于Java的部份,使用 NumberFormat來作格式化。