指定位數的圓周率

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

解法思路

首先介紹 J.Marchin 的圓周率公式：: PI = [16/5 - 16/(3*5³) + 16/(5*5⁵) - 16/(7*5⁷) + ……] - [4/239 - 4/(3*239³) + 4/(5*239⁵) - 4/(7*239⁷) + ……]
可以將這個公式整理為：: PI = [16/5 - 4/239] - [16/5³ - 4/239³]/3 + [16/5⁵ - 4/239⁵]/5 + ……
也就是說第 n 項，若 n 為奇數則為正數，為偶數則為負數，而項數表示方式為：: [16/5²ⁿ⁻¹ - 4/239²ⁿ⁻¹] / (2*n-1)
如果要計算圓周率至 10 的負 L 次方，由於 [16/5²ⁿ⁻¹ - 4/239²ⁿ⁻¹] 中 16/5²ⁿ⁻¹ 比 4/239²ⁿ⁻¹ 來得大，具有決定性，表示至少必須計算至第 n 項：: [16 /5²ⁿ⁻¹] / (2*n-1) = 10⁻ᶫ
將上面的等式取 log 後化簡，可以求得：: n = L/(2log5) = L/1.39794
若要求精確度至小數後 L 位數，則只要求至公式的第 n 項，其中 n 等於：: n = [L / 1.39794] + 1
在上式中 [] 為高斯符號，也就是取至整數（不大於 L/1.39794 的整數）；為了計算方便，可以在程式中使用下面這個公式來計算第 n 項：: [Wₙ₋₁/5²ⁿ - Vₙ₋₁/239²ⁿ] / (2*n-1)

其中 n 從 1 開始，而 W ₀ 為 16 * 5，V₀ 為 4 * 239，這個公式的演算法配合大數運算函式的演算法為：

div(w, 25, w);
divide(v, 57121, v); // 239 * 239 = 57121
sub(w, v, q);
div(q, 2 * k - 1, q);

最後的 PI 就是由單獨求得的第 n 項加總而得，因而大數運算後第一個數字代表整數，第二個位數之後是代表小數。為了計算精確度至小數後 L 位數，大整數運算必須有 L + 1 個位數。如果具備處理浮點數精確度的 API，亦可直接使用，但要注意預留位數與四捨五入問題，因為在除法過程中，可能有循環小數無法表示的問題。

程式實作：

#include <stdio.h> 
#include <stdlib.h>

#define L 1000 
#define N L / 4 + 1 

// L 為位數，N是array長度 

// 只處理正數的大整數加、減、除
void add(int*, int*, int*);       
void subtract(int*, int*, int*);
void divide(int*, int, int*);

int main(void) { 
    int s[N] = {0}; 
    int w[N] = {0}; 
    int v[N] = {0}; 
    int q[N] = {0}; 
    int n = (int) (L / 1.39793 + 1); 

    w[0] = 16 * 5; 
    v[0] = 4 * 239;

    int k; 
    for(k = 1; k <= n; k++) { 
        // 套用公式 
        divide(w, 25, w); 
        divide(v, 57121, v); // 239 * 239 = 57121
        subtract(w, v, q); 
        divide(q, 2 * k - 1, q); 

        if(k % 2) {// 奇數項 
            add(s, q, s); 
        } else {   // 偶數項 
            subtract(s, q, s); 
        }
    } 

    printf("%d.", s[0]); 
    for(k = 1; k < N; k++) {
        printf("%04d", s[k]); 
    }
   
    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 subtract(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 divide(int* a, int b, int *c) {  // b 為除數 
    int i, tmp, remain = 0; 
    for(i = 0; i < N; i++) { 
        tmp = a[i] + remain; 
        c[i] = tmp / b; 
        remain = (tmp % b) * 10000; 
    } 
}

import java.math.BigInteger;

public class LongPI {
    private final BigInteger PI;
    
    public LongPI(int L) {
        int n = (int) (L / 1.39793 + 1);
        
        BigInteger b25 = BigInteger.valueOf(25);
        BigInteger b57121 = BigInteger.valueOf(57121);
        BigInteger scale = BigInteger.valueOf(10).pow(L);
        
        BigInteger w = BigInteger.valueOf(16 * 5).multiply(scale);
        BigInteger v = BigInteger.valueOf(4 * 239).multiply(scale);
        BigInteger s = BigInteger.valueOf(0);
        BigInteger q = null;
        for(int k = 1; k <= n; k++) {
            w = w.divide(b25);
            v = v.divide(b57121);
            q = w.subtract(v).divide(BigInteger.valueOf(2 * k - 1));
            s = k % 2 == 1 ? s.add(q) : s.subtract(q);
        }
        PI = s;
    }
    
    public String toString() {
        String str = PI.toString();
        return String.format("%c.%s", str.charAt(0), str.substring(1));
    }
    
    public static void main(String[] args) {
        System.out.println(new LongPI(1000));
    }
}

class LongPI:
    def __init__(self, L):
        n = int(L / 1.39793 + 1)
        scale = 10 ** L
        
        def pi(k, w, v):
            if k == n + 1:
                return 0
            else:
                wk = w // 25
                vk = v // 57121
                qk = (wk - vk) // (2 * k - 1)
                return (qk if k % 2 == 1 else -qk) + pi(k + 1, wk, vk)
        
        self.PI = pi(1, 16 * 5 * scale, 4 * 239 * scale)
        
    def __str__(self):
        s = str(self.PI)
        return "%c.%s" % (s[0], s[1:])            
        
print(LongPI(1000))

import scala.math.BigInt

class LongPI private (pi: BigInt) {
    override def toString = {
        val str = pi.toString
        "%c.%s".format(str(0), str.tail)
    }
}

object LongPI {
    def apply(l: Int) = {
        val n = (l / 1.39793 + 1).toInt
        val b25 = BigInt(25)
        val b57121 = BigInt(57121)
        val scale = BigInt(10).pow(l)
        
        def pi(k: Int, w: BigInt, v: BigInt): BigInt = {
            if(k == n + 1) BigInt(0)
            else {
                val wk = w / b25
                val vk = v / b57121
                val qk = (wk - vk) / BigInt(2 * k - 1)
                (if(k % 2 == 1) qk else -qk) + pi(k + 1, wk, vk)
            }
        }
        
        new LongPI(pi(1, BigInt(16 * 5) * scale, BigInt(4 * 239) * scale))
    }
}

print(LongPI(1000))

class LongPI
    def initialize(l)
        n = (l / 1.39793 + 1).to_i
        scale = 10 ** l
        
        p = ->(k, w, v) {
            if k == n + 1
                0
            else
                wk = w / 25
                vk = v / 57121
                qk = (wk - vk) / (2 * k - 1)
                (if k % 2 == 1; qk else -qk end) + p.call(k + 1, wk, vk)
            end            
        }
        
        @pi = p.call(1, 16 * 5 * scale, 4 * 239 * scale)
    end
    
    def to_s
        str = @pi.to_s
        sprintf("%c.%s", str[0], str[1, str.size])
    end
end

puts(LongPI.new(1000))

// 用到大數運算中的 BigNumber

var LongPI = function() {
    function apply(L) {
        var n = parseInt(L / 1.39793 + 1);
        var b25 = BigNumber('25');
        var b57121 = BigNumber('57121');
        var scale = BigNumber('1' + new Array(L + 1).join('0'));
        
        var w = BigNumber('80').multiply(scale);
        var v = BigNumber('956').multiply(scale);
        var s = BigNumber('0');
        for(var k = 1; k <= n; k++) {
            w = w.divide(b25);
            v = v.divide(b57121);
            q = w.subtract(v).divide(BigNumber((2 * k - 1) + ''));
            s = k % 2 === 1 ? s.add(q) : s.subtract(q);
        }
        return new LongPI(s);
    }
    
    function LongPI(PI) {
        this.PI = PI;
    }
    
    LongPI.prototype.toString = function() {
        var str = this.PI.toString();
        return str[0] + '.' + str.substring(1);
    };
    
    return apply;
}();

print(LongPI(50));

data LongPI = LongPI Integer

instance Show LongPI where
    show (LongPI value) = (head s) : '.' : (tail s)
        where s = show value

longPi l = LongPI $ p 1 (16 * 5 * scale) (4 * 239 * scale)
    where n = truncate (fromIntegral l / 1.39793 + 1)
          scale = 10 ^ l
          p k w v =
              if k == n + 1 then 0
              else
                  let wk = w `div` 25
                      vk = v `div` 57121
                      qk = (wk - vk) `div` (2 * k - 1)
                  in (if odd k then qk else -qk) + (p (k + 1) wk vk)

main = print $ longPi 1000