多倍長精度による円周率 π の計算

「多倍長精度による計算」に，R による円周率 π の計算プログラムが掲載されている。R でも gmp パッケージを使えば多倍長精度による計算はできるので，このプログラムを書き換えればよい。

このプログラムは円周率を10の羃乗倍した固定小数点数として計算する。

（浮動小数点数で計算するプログラムは「多倍長精度計算パッケージ Rmpfr」に掲載した。）

末尾に Python と Julia のプログラムを追記した。

library(gmp)
calc.pi = function(accuracy=10000)    {
  FACTOR = as.bigz(10)^accuracy # FACTOR 倍した固定小数点数として計算
  pi = as.bigz(0)
  a = as.bigz(16*5*FACTOR)
  n = 1
  repeat {
    a = a %/% 25
    b = a %/% n
    if (b == 0) break
    pi = pi + b
    n = n + 2
    a = a %/% 25
    b = a %/% n
    if (b == 0) break
    pi = pi - b
    n = n + 2
  }
  a = as.bigz(4*239*FACTOR)
  n = 1
  repeat {
    a = a %/% 57121
    b = a %/% n
    if (b == 0) break
    pi = pi - b
    n = n + 2
    a = a %/% 57121
    b = a %/% n
    if (b == 0) break
    pi = pi + b
    n = n + 2
  }
  return(pi)
}

> calc.pi(50)
Big Integer ('bigz') :
[1] 314159265358979323846264338327950288419716939937510

Python でも gmpy モジュールを使えば多倍長精度による計算ができる。

移植に際して注意すべきことは，除算は / ではなく，整数除算 // を使うこと。また，π は浮動小数点数であるが，指定した桁の固定小数点数として扱う。具体的には小数点以下の桁数（acculacy で指定）の下に小数点を置いて，整数値 π x acculacy を計算する。

def calc_pi(acculacy=10000):
    FACTOR = 10**acculacy
    pi = 0
    n, a = 1, 16 * 5 * FACTOR
    while True:
        a //= 25
        b = a // n
        if b == 0:
            break
        pi += b
        n += 2
        a //= 25
        b = a // n
        if b == 0:
            break
        pi -= b
        n += 2
    n, a = 1, 4 * 239 * FACTOR
    while True:
        a //= 57121
        b = a // n
        if b == 0:
            break
        pi -= b
        n += 2
        a //= 57121
        b = a // n
        if b == 0:
            break
        pi += b
        n += 2
    return pi

>>> calc_pi(50)
314159265358979323846264338327950288419716939937510

実際の小数点は，言うまでもなく，先頭桁の 3 の次にある。

なお，Python には，sympy があり，それを用いれば，任意の精度の π が簡単に得られる。

>>> import sympy
>>> sympy.pi.evalf(50)
3.1415926535897932384626433832795028841971693993751

50 桁までの指定では，sympy.pi と比べると，最後まで合っている。普通は，最後の2桁ぐらいに誤差が出る。

Julia の場合も裏で gmp ライブラリを呼んでいるので，ほぼ同じように書くことができる。
多倍長精度の整数は big(数) で定義し，整数除算は ÷，if (b == 0) break が b == 0 && break と書けるなど。

function calcpi(accuracy = 10000)
    factor = big(10) ^ accuracy
    pi = big(0)
    a = 16 * 5 * factor
    n = big(1)
    while true 
        a ÷= 25
        b = a ÷ n
        b == 0 && break
        pi += b
        n += 2
        a ÷= 25
        b = a ÷ n
        b == 0 && break
        pi -= b
        n += 2
    end
    a = 4 * 239 * factor
    n = big(1)
    while true 
        a ÷= 57121
        b = a ÷ n
        b == 0 && break
        pi -= b
        n += 2
        a ÷= 57121
        b = a ÷ n
        b == 0 && break
        pi += b
        n += 2
    end
    pi
end