From Gossip@caterpillar

Algorithm Gossip: 長 PI

說明

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

解法

首先介紹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來作格式化。

您可以參考 這個網站 有關於另一個用公式求PI的說明，它的實作在 這邊。

實作

  • C 
#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; 
    } 
}

  • Java 
import java.text.NumberFormat;

public class LongPI {
    public static void main(String[] args) {
        final int L = 1000;
        final int N = L/4+1;
        int[] s = new int[N+3];
        int[] w = new int[N+3];
        int[] v = new int[N+3];
        int[] q = new int[N+3];

        int n = (int)(L/1.39793 + 1); 
        
        w[0] = 16*5; 
        v[0] = 4*239; 
        
        for(int k = 1; k <= n; k++) { 
            // 套用公式 
            w = BigNumber.div(w, 25); 
            v = BigNumber.div(v, 239); 
            v = BigNumber.div(v, 239);
            q = BigNumber.sub(w, v); 
            q = BigNumber.div(q, 2*k-1); 

            if(k % 2 == 1) // 奇數項 
                s = BigNumber.add(s, q); 
            else    // 偶數項 
                s = BigNumber.sub(s, q);
        } 
        
        NumberFormat nf = NumberFormat.getInstance();
        nf.setMinimumIntegerDigits(4);
        
        System.out.print(s[0] + "."); 
        for(int k = 1; k < N; k++) 
            System.out.print(nf.format(s[k])); 
        
        System.out.println(); 
    }
}