💫

任意精度演算を使って円周率を計算する

2024/06/23に公開

Julia

tech

背景

新しい円周率の級数表現が流行していたので、Juliaの任意精度演算を使って計算してみました。

基本的な使い方

基本的な機能はサポートされているので、setprecisionくんを使います。64bitと256bitで定数 $\pi$ を出力してみた例です。だいたい長さが4倍ぐらいになっていることが見てとれます。

# BigFloat演算のビット数 (= 64bit)
julia> setprecision(64) do
           println(BigFloat(π))
       end
3.14159265358979323851

# 256bitに設定
julia> setprecision(256) do
           println(BigFloat(π))
       end
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

公式のドキュメントも参照ください（任意精度演算）

実装

先ほど確認した任意精度演算の機能を使って、若井さんのPython実装をJuliaで動かしてみます。

基本的な計算

パラメータZで精度を指定してから級数を計算します。

Pythonのmpmathのnsumでは級数計算が高性能に実装されている様子でしたが、適当にmaxiterを指定して計算しました。
複素数にも対応できるように出力時にreal()を取りましたが、今回はこの機能は使いませんでした。
Pochhammerのべき乗計算のところで、0とn-1あたりを最初実装ミスして $\pi$ に漸近せず、時間がかかりました（1敗）

function compute_pi(λ;Z=512, maxiter=1024)
    setprecision(Z)

    function Pochhammer(x, n)
        result = big(1)
        for i in 0:n-1
            result *= x + big(i)
        end
        result
    end

    function f(n, λ)
        x = ((2n + 1) ^ 2 / (4 * (n + λ))) - n
        (1 / (n + λ) - 4 / (2n + 1)) * Pochhammer(x, n - 1) / factorial(n)
    end

    real(4 + sum(f(big(k), λ) for k in 1:maxiter))
end

計算した数値が、Juliaの持っている $\pi$ と比較し、どの桁が一致するかをprintする関数も書きました。

using Printf
function print_diff(λ;Z=512, maxiter=1024)
    setprecision(Z)
    str_cpi = "$(compute_pi(λ;Z=Z, maxiter=maxiter))"
    str_pi = "$(BigFloat(π))"

    println("   j | 1 2")
    println("-----------")
    for j in 1:length(min(str_cpi, str_pi))
        flag = str_cpi[j] != str_pi[j] ? "*" : ""
        head = @sprintf("%4d", j)
        println("$(head) | $(str_cpi[j]) $(str_pi[j]) $flag")
    end
end

実装したものを実行してみました ( $\lambda=10$ )。実行結果を見ると、31桁目までは一致していたようでした。

julia> print_diff(10)
   j | 1 2
-----------
   1 | 3 3 
   2 | . . 
   3 | 1 1 
   4 | 4 4 
   5 | 1 1 
   6 | 5 5 
   7 | 9 9 
   8 | 2 2 
   9 | 6 6 
  10 | 5 5 
  11 | 3 3 
  12 | 5 5 
  13 | 8 8 
  14 | 9 9 
  15 | 7 7 
  16 | 9 9 
  17 | 3 3 
  18 | 2 2 
  19 | 3 3 
  20 | 8 8 
  21 | 4 4 
  22 | 6 6 
  23 | 2 2 
  24 | 6 6 
  25 | 4 4 
  26 | 3 3 
  27 | 3 3 
  28 | 8 8 
  29 | 3 3 
  30 | 2 2 
  31 | 7 7 
  32 | 6 9 *

誤差の変化

$\pi$ との誤差がどのように変化するのかも確認しました。次のような実験用の関数を作りました。

$\lambda$ を10から300まで10刻みで変化させる
精度 $Z$ を受取指定する
最小の反復回数 maxiter と倍率 iterK を受取、maxiter、2*maxiter、… と変化させた結果を計算する
Juliaの持っている $\pi$ との誤差 (の絶対値) を対数プロットする
- ここでたまに0が出てしまってグラフ描くのに不都合があったので、適当に1e-128足してます

function diff_maxiter(;Z=2048, maxiter=128, iterK=6)
    R = 10:10:300
    
    f = plot(size=(500, 300), dpi=150)

    for k in 1:iterK
        mik = maxiter * k
        diff = setprecision(Z) do
            K = zeros(BigFloat, length(R))
            for (ir, λ) in enumerate(R)
                pi_λ = compute_pi(λ; Z=Z, maxiter=mik)
                K[ir] = pi_λ
            end
            abs.(BigFloat(π) .- K) .+ BigFloat(1e-128)
        end
        plot!(f, R, diff, label="loss(Z=$Z, mi=$mik)", yscale=:log10)
    end

    savefig(f, "figures/loss_diff_mi_Z$(Z)_K$(iterK).png")
end