# 長 PI

## 解法

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 + ......

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

[16 / 5 2 * n - 1 ] / (2 * n - 1) = 10-L

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

n = [L / 1.39794] + 1

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

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

## 實作：CJavaPythonScalaRubyJavaScriptHaskell

• C
``#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.BigIntclass 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])    endendputs(LongPI.new(1000))``

• JavaScript
``// 用到 大數運算 中的BigNumbervar 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 Integerinstance Show LongPI where    show (LongPI value) = (head s) : '.' : (tail s)        where s = show valuelongPi 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``