「多倍長精度による計算」に，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
