長 PI

說明

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

解法

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

[16 / 5^{2 * n - 1} - 4 / 239^{2 * n - 1}] / (2 * n - 1)

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

[16 / 5^{2 * 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項：

[W_{n - 1} / 5²ⁿ - V_{n - 1} / 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個位數。如果具備處理浮點數精確度的BigDecimal之類的API，亦可直接使用，但要注意預留位數與四捨五入問題，因為在除法過程中，可能有循環小數無法表示的問題。

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

實作：C Java Python Scala Ruby JavaScript Haskell

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

Java

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));
    }
}

Python

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))

Scala

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))

Ruby

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))

JavaScript

// 用到 大數運算 中的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));

Haskell

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