Juliaで円周率10万桁を解析して数字の規則性を可視化してみた

はじめに

新しいプログラミング言語に挑戦したくて、Julia言語でのコーディング学習を始めました。

Juliaは、任意精度の浮動小数点数 BigFloat が標準で使えたり、グラフ描画のパッケージが充実していたりと、数値実験にぴったりの環境です。

せっかくなので、円周率πの小数部分で現れる数字に偏りはあるのか について調べてみることにしました。

やったこと

• πの小数点以下の各数字（0〜9）の出現頻度を集計
• 連続する2桁ペア（00〜99）の出現頻度も集計
• 棒グラフ・折れ線グラフ・ヒートマップで可視化

環境

• Julia 1.6
• 使用パッケージ: Plots

πの小数部分を取得する

任意精度の設定

setprecision(BigFloat, 332192)  # 約10万桁
pi_big = BigFloat(π)

setprecision(BigFloat, n) は、BigFloat の精度を ビット単位 で指定する関数です。BigFloat(π) と書くだけで指定した精度のπが得られるのは、凄いですよね。

なぜ332,192なのか

最初は setprecision(BigFloat, 100000) と書いて10万桁を期待していたのですが、実際に取得できた小数部分の長さは 約30,103桁 でした。「10万桁のつもりだったのに……」と戸惑いましたが、これは setprecision の引数が ビット数 であり、十進の桁数ではない ためです。

ビット数から十進桁数への変換は、以下の式で求められます。

$$
\text{十進桁数} \approx \frac{\text{ビット数}}{\log_2 10}
$$

$\log_2 10 \approx 3.3219$ なので、

$$
\frac{100000}{3.3219} \approx 30103
$$

つまり、100,000ビットの精度は十進で約30,103桁に相当する、ということです。

10万桁を得るには、逆算して必要なビット数を求めます。

$$
\text{必要ビット数} = \lceil 100000 \times \log_2 10 \rceil = \lceil 100000 \times 3.3219... \rceil = 332193
$$

332193を設定したら少しずれたので、調整し332192 ビットを指定しました。

setprecision(BigFloat, 332192)  # これで約10万桁が得られる

小数部分の文字列化

pi_str = string(pi_big)[3:end]  # "3." の2文字を除去

string(pi_big) は "3.14159265..." のような文字列を返すので、先頭の "3." をスライスで取り除いて小数部分だけにしています。Juliaの文字列インデックスは1始まりなので、[3:end] で3文字目以降を取得します。

各数字の出現回数をカウント

counts = Dict(digit => 0 for digit in '0':'9')
for c in pi_str
    counts[c] += 1
end

Dict の内包表記で '0' 〜 '9' をキーとする辞書を作り、1文字ずつループしてカウントしています。

出力結果

小数部分の長さ: 100000
数字の出現回数:
0: 9999
1: 10137
2: 9908
3: 10025
4: 9971
5: 10026
6: 10028
7: 10026
8: 9978
9: 9902

可視化

棒グラフ

observed = [counts[c] for c in '0':'9']
bar(string.('0':'9'), observed,
    title="Digit Distribution of π (100,000 digits)",
    xlabel="Digit", ylabel="Count",
    legend=false)
savefig("bar_chart.png")

縦軸のスケールに注目すると、ほぼ均等に分布していることがわかります。

累積カウントの折れ線グラフ

桁数が増えるにつれて各数字のカウントがどう推移するかを可視化します。もし偏りがあれば、特定の線だけが大きく離れるはずです。

n = length(pi_str)
cumcounts = zeros(Int, n, 10)  # n行×10列の行列
running = zeros(Int, 10)       # 各数字の現在の累積カウント
for i in 1:n
    d = pi_str[i] - '0'  # 文字を整数（0〜9）に変換
    running[d+1] += 1     # Juliaは1始まりなので+1
    cumcounts[i, :] .= running  # i行目に現在の累積値をコピー
end

plot(1:n, cumcounts,
    title="Cumulative Digit Count of π",
    xlabel="Number of digits",
    ylabel="Count",
    label=reshape(string.(0:9), 1, 10),
    legend=:topleft,
    size=(800, 500))
savefig("line_chart.png")

10本の線がほぼ重なり合いながら直線的に伸びていて、特定の数字だけが突出しているという傾向は見られません。

2桁ペアの解析

1桁だけでなく、連続する2桁（00〜99）についても同じ分析を行ってみました。

スライディングウィンドウでペアを取得

pairs = [pi_str[i:i+1] for i in 1:length(pi_str)-1]

例えば小数部分が "14159" なら、["14", "41", "15", "59"] というペアが得られます。隣接する2文字を1文字ずつずらしながら抽出する、いわゆるスライディングウィンドウ方式です。

カウント

pair_counts = Dict(lpad(i, 2, '0') => 0 for i in 0:99)
for p in pairs
    pair_counts[p] += 1
end

pair_labels = [lpad(i, 2, '0') for i in 0:99]
pair_observed = [pair_counts[k] for k in pair_labels]

lpad(i, 2, '0') は数値 i を2桁にゼロ埋めする関数です（例: 7 → "07"）。

2桁ペアの出現回数

2桁ペアの出現回数:
00: 998  01: 1027  02: 962  03: 993  04: 968  05: 1007  06: 1009  07: 1017  08: 1001  09: 1017
10: 1042  11: 1034  12: 992  13: 1010  14: 1030  15: 971  16: 982  17: 1008  18: 1044  19: 1024
20: 974  21: 1064  22: 971  23: 966  24: 943  25: 1003  26: 1048  27: 1017  28: 955  29: 967
30: 979  31: 974  32: 994  33: 1008  34: 1041  35: 1075  36: 965  37: 1009  38: 982  39: 998
40: 1020  41: 1001  42: 987  43: 1014  44: 971  45: 962  46: 1032  47: 1009  48: 975  49: 1000
50: 966  51: 966  52: 1029  53: 1031  54: 1019  55: 1015  56: 1017  57: 1004  58: 1008  59: 971
60: 1054  61: 1012  62: 1048  63: 997  64: 1020  65: 998  66: 953  67: 972  68: 949  69: 1025
70: 948  71: 1046  72: 981  73: 1031  74: 945  75: 961  76: 1100  77: 1012  78: 1040  79: 961
80: 1011  81: 996  82: 992  83: 946  84: 1001  85: 1020  86: 970  87: 1044  88: 1027  89: 971
90: 1007  91: 1016  92: 952  93: 1029  94: 1033  95: 1014  96: 952  97: 934  98: 997  99: 968

2桁ペアの棒グラフ

bar(pair_labels, pair_observed,
    title="2-Digit Pair Distribution of π",
    xlabel="Pair", ylabel="Count",
    legend=false,
    xticks=(1:5:100, pair_labels[1:5:100]),
    xrotation=45,
    size=(1000, 500))
savefig("pair_bar_chart.png")

法則性があるかもしれませんが、少しばらつきがあるくらいで100種類のペアがほぼ均一に出現しています。

ヒートマップ

heatmap_data = reshape(pair_observed, 10, 10)'  # 行:十の位, 列:一の位
heatmap(string.(0:9), string.(0:9), heatmap_data,
    title="2-Digit Pair Frequency of π",
    xlabel="Second digit",
    ylabel="First digit",
    color=:viridis,
    size=(600, 500))
savefig("pair_heatmap.png")

reshape(pair_observed, 10, 10) で、100要素の1次元配列を10×10の行列に変形しています。行が十の位、列が一の位に対応します。

全体的に均一な色合いで、特定のペアが極端に多い・少ないという傾向は見られません。

まとめ

今回はJulia言語を用いて、円周率の小数部分に規則性があるかを可視化してみました。

結果としては、規則性は見つからない ということが分かりました！！

もしコンピューターがなければ、円周率を手計算し、その数字をひとつひとつ目で追ってカウントしていく作業になります。数か月かけた挙句に規則性が見つからないなんてことも十分あり得ます。それをわずか数時間で検証できるのは、プログラミングの力は偉大ですね。

また、Julia言語の特徴である計算速度も、10万桁の解析が実行から約10秒でグラフ作成まで完了し、速さを体感できました。引き続き、勉強していきたいと思います。

下記が今回のGithubのリポジトリになります。興味のある方はどうぞ！

ここまで読んでいただきありがとうございました。