あなたは円周率を何桁言えますか。3.14159…という、あの数字です。 円周率の小数部分は無限に続き、循環することもありません。 古来より、数学者は円周率の値を様々な幾何学的な近似や公式を用いて計算してきました。 その桁数は計算機の発明により飛躍的に伸び、収束の速い公式の発見や効率の良いアルゴリズムの発明などによって加速してきました *1。
5年前、私がまだ学生だった頃、円周率1億桁の計算に挑んだことがありました。 私にとって高精度計算の初めての挑戦で、様々な試行錯誤で苦労したのをよく覚えています。
2017年現在、円周率計算の世界記録は22兆桁です。 円周率計算の歴史をご覧いただくとよく分かると思いますが、近年の円周率計算の世界記録からは次のような特徴が読み取れます。
- 2002年に1兆を超え、最新の記録 (2016年) は22兆桁 (10進数)
- 企業ではなく個人によって記録が塗り替えられてきている
- Alexander J. Yee氏によって作られたy-cruncherが多くの記録で使われている
- ストレージは100TB単位
- 計算時間は100日〜1年がかり
追記: 2019年3月14日、Googleによって31兆桁 (31415926535897桁) に記録が更新されました。やはり同じくy-cruncherを使っており、検証アルゴリズムもBellard's formulaとBBP formulaの両方を使ったとのことです (参考)。自社のクラウドリソースを自由に使えるビッグプレイヤーが記録を更新してきた以上、今後個人でこれに立ち向かうのは厳しいかと思われます。
世界記録が22兆桁である一方で、5年前に私が計算したのはたったの1億桁です。 当時、私のブログのコードを少し改変して10億桁計算された方もいらっしゃいました (ありがとうございます) が、それでも世界記録とは1万倍の差があります。 そこには、かけた計算時間以上の技術の差があります。
- 数字の桁数が多すぎてメモリーに乗らないので、ファイルに書きつつ計算する必要がある
- 掛け算1つするにしても、桁を分割してファイルを読み書きしながら行う必要がある
- そういうファイルに書きつつ計算するソフトウェアを自前で実装したり工夫する必要がある (億桁はGMPで計算できてしまった)
- 世界記録の桁は他の人が計算したことがない値なので、計算結果を検証する方法を工夫する必要がある (私はオンラインサービスの結果と突き合わせてしまった)
本記事では、4番目の、円周率計算の結果の検証について説明してみようと思います。 円周率計算の記録は、どうやって検証するのでしょうか。 渋いテーマですが、世界記録をどう検証するかには興味がありますし、5年前の自分がスルーしてしまった部分でした。 当時の私は、オンラインで任意の桁を表示できるページで答えを見て済ませてしまっていました。 それができたのは、1億桁程度だからだったのと、たまたまそういうオンラインサービスがあったからです。 「人類未踏の値」の検証の裏には、数学と計算機科学の織りなす美しい公式の発見がありました。
- 作者:牧野 貴樹
- 発売日: 1996/03/01
- メディア: 単行本
円周率の世界記録の検証
円周率の世界記録が正しいことを検証するにはどうすればいいのでしょうか。 求めた桁の値は、その計算をした人しか見たことがない値です。 別の公式を使って同じ桁数計算して、結果を突き合わせるという方法がありますが、これには膨大な時間がかかります。 20兆桁計算するのに100日ほどかかっているので、別のアルゴリズムで全く同じ桁数分計算するにはさらに100日以上かかってしまいます。
ところが円周率計算の歴史を見ていると、最新の世界記録の検証は28時間で済んでいると書かれています。 2014年の記録でも182時間、その前年の記録でも46時間となっています。 全体の掛かった時間が100日という単位ですから、それに比べるとこれらの検証時間は短すぎます。 この時間で、別のアルゴリズムで全ての桁数計算したとは思えません。 もしこの時間で全く同じ桁数の計算を行えるなら、最初からそれを使ってますよね。 これは一体どういうことなのでしょうか。
実は、全ての桁数の検証をしているわけではないのです。 円周率の公式でBBP formulaと呼ばれているものがあります *2。
この公式を使えば、円周率の16進数表現の任意の桁を、高速に求めることができるのです。
つまり16進数表現の小数点以下 桁目の数を、
桁までの数字を使うことも保存することもなく、直接求めることができるのです。
円周率計算の世界において、この公式が発見されたことは衝撃的でした。
今から当時の衝撃を感じるのは困難なのですが、それ以前はこの公式がなく、愚直に別のアルゴリズムで検証する他なかったことを想像すると、とても重要な公式の発見だと言えるでしょう。
BBP formulaの肝は、 の分母が2の累乗であることです。
つまり、この公式を使うと、計算機上の数字の表現の生のデータを、直接検証することができます。
実はBBP公式の発見以降、類似の公式がいくつか見つかっています。
そのうちの1つが1997年に発見されたBellard's formulaです *3。
この公式も分母が2の累乗なので、円周率の16進数表現の検証に用いることができます。 同じ目的の公式が2つ以上あることで、仮に一方の計算に誤りがあったとしても、もう片方でさらに検証を重ねることができます。
以上のことを踏まえて、円周率計算の世界記録がどう計算され、どう検証されるのかをまとめてみました。
- 円周率の計算を行う。2進数表現のままファイルに書いたり読んだりしながら計算が進みます。兆桁台にもなると、計算中にハードディスクが壊れるといったリスクもありますので、復帰ポイントを随時作りながら、プログラムをいったん止めても再開できるようにしておく必要があります。近年の世界記録はChudnovskyアルゴリズム*4にもとづいています。1989年に発見されたChudnovsky兄弟によって発見された公式で、収束が高速である (一項計算すると14桁進む) ことが知られています。
- 2進数表現の円周率を、10進数表現に変換します。
- 2進数表現の検証を行います。円周率の最後の方の桁を、BBP formulaやBellard's formulaのような任意桁を求める公式によって検証します。イテレーションの早い段階や浅い桁で起きたエラーはすぐに最後の桁まで伝搬するので、この検証でかなり正当性が保証されます。
- しかし前の桁で間違っている可能性は否定できません。特に、最後の掛け算中に浅い桁で誤りが起きる可能性があります。最終フェーズの掛け算は剰余を使って検証します。つまり十分大きな素数
に対して、
を確認します。
- 10進数表現の検証を行います。2進数表現が正しくても基数変換に誤りがある可能性があるので、ここも検証する必要があります。10進数表現を2進数表現に変換して元に戻れば確実です。素数で割った余りを比較する方法もあります。つまり素数
に対して
を計算する方法です。10進数は10進数のまま、16進数は16進数のまま10掛けたり剰余を計算し、最後に計算を突き合わせるのです。
整数演算は2進数で行いますが、最後は10進数表現の円周率に変換しないと何兆桁計算したとはいえません。 この基数変換の検証も重要な部分なのですが、このエントリーでは省略します。 なお、10進数版BBP formulaがもし発見できれば、円周率10進数表現の任意桁をすぐに求められるようになり、計算結果の検証どころか、何百兆桁, 何千兆桁目を現実的な時間の中で求めることさえできてしまうという驚異的な結果なのですが、2017年現在ではそういう式は見つかっていません。
円周率計算の検証の更なる詳細はPi - 5 Trillion DigitsやAlexander Yeeさんのstackoverflowの回答、2016年世界記録のブログなどをご覧ください。 BBP formulaやBellard's formulaがいかに重要かがよく分かると思います。
BBP formula
BBP formulaを実装してみましょう。 あらためて公式を眺めてみます。
求めたいのは円周率の16進数表現の、小数点以下 桁目の数字です。
ここで小数部分を返す
(
はガウス記号で
以下の最大の整数) を定義すると、
と書くことができます。ただし
と定義しました。 ここで和を2つに分けて
と書けますが、一項目の は整数で、かつ小数部分を返すのは和のそれぞれの項に分けても同じ (
) を用いると次のように書くことができます。
第一項の は何千万何億となるので、
が大きくなりすぎないように剰余を取るのが大事なのです (いわゆる「冪剰余」です)。
第二項はどんどん0に収束していくので、適当なところまで計算したら和を取ってしまいます。
近年Rustが大流行していて、利用者が爆発的に増加しているという状況の中で、私も早く入門したいという気持ちが強まっています。 円周率計算はちょうどよい題材なのでRustで実装してみました。
単一ファイルで実行できるように、ここにもおいておきます。
pihex_bbp.rs
use std::env; use std::thread; use std::sync::mpsc; fn main() { let place: u64 = env::args().nth(1).and_then(|x| x.parse().ok()).unwrap_or(0); print!("{}:", place); for i in 0..4 { print!(" {}", pihex_bbp(place + 4 * i)); } print!("\n"); } fn pihex_bbp(d: u64) -> String { let (tx, rx) = mpsc::channel(); for &(j, k) in &[(1, 4.0), (4, -2.0), (5, -1.0), (6, -1.0)] { let tx = tx.clone(); thread::spawn(move || tx.send(k * series_sum(d, j)).unwrap()); } drop(tx); let fraction: f64 = rx.iter().sum(); (0..4) .scan(fraction, |x, _| { *x = (*x - x.floor()) * 16.0; Some(format!("{:x}", x.floor() as u64)) }) .fold(String::new(), |s, t| s + &t) } fn series_sum(d: u64, j: u64) -> f64 { let fraction1: f64 = (0..d + 1) .map(|i| (pow_mod(16, d - i, 8 * i + j) as f64) / ((8 * i + j) as f64)) .fold(0.0, |x, y| (x + y).fract()); let fraction2: f64 = (d + 1..) .map(|i| 16.0_f64.powi(-((i - d) as i32)) / ((8 * i + j) as f64)) .take_while(|&x| x > 1e-13_f64) .sum(); fraction1 + fraction2 } fn pow_mod(n: u64, m: u64, d: u64) -> u64 { if m == 0 { 1 % d } else if m == 1 { n % d } else { let k = pow_mod(n, m / 2, d); if m % 2 == 0 { (k * k) % d } else { (k * k * n) % d } } }
実行してみましょう
$ rustc -O pihex_bbp.rs $ ./pihex_bbp 1000000 1000000: 6c65 e52c b459 3500 # ./pihex_bbp 1000000 4.52s user 0.04s system 305% cpu 1.494 total $ ./pihex_bbp 10000000 10000000: 7af5 863e fed8 de97 # ./pihex_bbp 10000000 52.90s user 0.30s system 351% cpu 15.154 total $ ./pihex_bbp 100000000 100000000: cb84 0e21 926e c5ae # ./pihex_bbp 100000000 606.44s user 3.26s system 361% cpu 2:48.52 total $ sysctl -n machdep.cpu.brand_string Intel(R) Core(TM) i7-5557U CPU @ 3.10GHz
およそ3分で1億桁目の数字が表示されました。
このコードでは4回pihex_bbpの計算を行っているため、4倍時間がかかっていることに注意してください。
1億桁目付近ののcb84は45秒ほどで計算することができます。
円周率自体の計算は、桁が増えるほど多くのメモリーやディスクを消費します。 一方で、BBP formulaを使ったアルゴリズムでは計算過程で桁が増えていくことはありません。 小数部分をdoubleでひたすら計算していくので、消費メモリーが増えていくこともありませんし、桁数が増えるために基本演算の計算時間が伸びていくといったこともありません。 これが、100日かかる世界記録をたった数十時間で検証できる理由なのです (実際はこの公式以外にも複合的な観点から検証されます。いずれにせよBBP formulaを使ってランダムサンプリングした桁の数字をチェックしたり、剰余を取って検証したりという感じで、計算結果の数字をくまなくチェックすることはありません)。
Bellard's formula
BBP formulaと並んでよく使われる公式がBellard's formulaです。 BBP formulaを使って検証したとしても、BBPの実装自体が誤っている可能性も否定できませんので、Bellard's formulaも実装して結果を突き合わせておきましょう。
公式はこのようになっています。
ここで
とおくと、次の関係がなりたちます。
この式を使って実装してみました。
pihex_bellard.rs
use std::env; use std::thread; use std::sync::mpsc; fn main() { let place: u64 = env::args().nth(1).and_then(|x| x.parse().ok()).unwrap_or(0); print!("{}:", place); for i in 0..4 { print!(" {}", pihex_bellard(place + 4 * i)); } print!("\n"); } fn pihex_bellard(d: u64) -> String { let (tx, rx) = mpsc::channel(); for &(j, k, l) in &[(4, 1, -32.0), (4, 3, -1.0), (10, 1, 256.0), (10, 3, -64.0), (10, 5, -4.0), (10, 7, -4.0), (10, 9, 1.0)] { let tx = tx.clone(); thread::spawn(move || tx.send(l * series_sum(d, j, k)).unwrap()); } drop(tx); let fraction: f64 = rx.iter().sum(); (0..4) .scan(fraction, |x, _| { *x = (*x - x.floor()) * 16.0; Some(format!("{:x}", x.floor() as u64)) }) .fold(String::new(), |s, t| s + &t) } fn series_sum(d: u64, j: u64, k: u64) -> f64 { let fraction1: f64 = (0..(2 * d + 2) / 5) .map(|i| { ((if i % 2 == 0 { 1.0 } else { -1.0 }) * pow_mod(4, 2 * d - 3 - 5 * i, j * i + k) as f64) / ((j * i + k) as f64) }) .fold(0.0, |x, y| (x + y).fract()); let fraction2: f64 = ((2 * d + 2) / 5..) .map(|i| -(-4.0_f64).powi(-((5 * i + 3 - 2 * d) as i32)) / ((j * i + k) as f64)) .take_while(|&x| x.abs() > 1e-13_f64) .sum(); fraction1 + fraction2 } fn pow_mod(n: u64, m: u64, d: u64) -> u64 { if m == 0 { 1 % d } else if m == 1 { n % d } else { let k = pow_mod(n, m / 2, d); if m % 2 == 0 { (k * k) % d } else { (k * k * n) % d } } }
実行してみます。
$ rustc -O pihex_bellard.rs $ ./pihex_bellard 1000000 1000000: 6c65 e52c b459 3500 # ./pihex_bellard 1000000 3.40s user 0.02s system 341% cpu 1.003 total $ ./pihex_bellard 10000000 10000000: 7af5 863e fed8 de97 # ./pihex_bellard 10000000 39.30s user 0.23s system 362% cpu 10.905 total $ ./pihex_bellard 100000000 100000000: cb84 0e21 926e c5ae # ./pihex_bellard 100000000 447.77s user 2.75s system 355% cpu 2:06.71 total
先程のBBPの実行結果と見比べると、およそ25%速くなっていることがわかります。 Bellard's formulaはBBP formulaよりも43%速いと言われます *5。 そこまでの記録は出ませんでしたが、確かにBBPより速いようだということは観測できました。
1000万桁目の悲劇
上記のコードは最初からこの形をしていたわけではありませんでした。 最初修正に手間取った部分は以下のコードでした。
let fraction1: f64 = (0..d + 1) .map(|i| (pow_mod(16, d - i, 8 * i + j) as f64) / ((8 * i + j) as f64)) .fold(0.0, |x, y| (x + y).fract());
この最後の行は元々次のように書いていました。
.sum();
当初、総和を取っているから.sum()でという短絡的な考えでコードを書いていました。
その結果1000万桁目は次のようになっていました。
$ pihex 10000000 10000000: 7af5 863f fed8 de97
この結果が間違っているのがわかりますか？ 次の結果が正解です。
10000000: 7af5 863e fed8 de97
sumを使うと、863eが863fとなってしまいます。
偶然The BBP Algorithm for Pi (David H. Bailey) に計算結果が掲載されており、それと見比べていた時に気がついたのです。
863efeという数字の並びからお察しの通り、誤差が積もってこういう間違いが起きてしまっています。
計算に用いた式をあらためて見てみましょう。
この第一項の総和に注目してみます。 各項は剰余を取るために0以上1以下の数なのですが、何千万項も足しているために、結果が何百万という数になって、それだけで6桁も損してしまいます。 総和を取ってから小数部分を返すのではなく、各項を足しては小数部分にすることを繰り返すことで、常に和を1以下に保ち、小数点以下の精度が落ちることを防ぐ必要があったのです。
エラーが起きやすい項が1000万桁付近にあったことはとても幸運でした。 計算結果が載っている論文があったことも恵まれていました。 もし仮にこの誤差に気がついていなければ、誤った実装のまま満足してしまったことでしょう。
総和を一気にとるのではなく、足しては小数部分にすることを返すことでかなり誤差は減らすことができました。
しかし厳密に言うと、これでもあらゆる桁に正しく返せるとは言えません。
倍精度浮動小数点の仮数部、つまり52bit以上0やfが揃っていたら、どうしても浮動小数点計算誤差の影響を受けてしまいます。
こういう場所が円周率のどの桁にあるかは知りませんが、そういう場所がないとは言いきれない以上、浮動小数点を使う計算で以って絶対に正しく答えを返すという保証はないのです。
BBP formulaやBellard's formulaを計算機で実装すると、浮動小数点の誤差を必ず受けうることになります。
上記のように小数部分のみで計算進めるのは必要な対処ですが、その対処をしても避けられないケースは存在するのです。
しかし、このアルゴリズムを使う意図を思い出してみましょう。
計算した円周率の最後の方の桁の検証のために、特定の桁を計算するために使います。
計算結果を見て0やfが続いているかどうか、すなわち誤差の影響を受けうるかどうかは目で見て分かるので、そこまで大きな問題にはならないのです。
10億桁目の悪夢
上記の浮動小数点計算誤差に起因するバグを修正したことで、かなりのエラーは減りました。 ところが、いざ10億桁目を計算してみると、まためちゃくちゃな答えを返すようになりました。 これくらいの桁にもなると計算自体に数分かかるようになり、デバッグも困難を極めました。
上記のコードで10億桁目を計算した結果、次のようになっていました (時間がかかるのでpihex_{bbp,bellard}を1回にしています)。
$ ./pihex_bbp 1000000000 1000000000: 4fa1de2c22c1 $ ./pihex_bellard 1000000000 1000000000: a3222cec6720
両者の結果が異なっている時点で少なくとも一方が間違っていると分かりますし、そもそも正解は5895585a0428...ですから (David H. Bailey, The BBP Algorithm for Pi)、両者ともめちゃくちゃです。
BBP formulaとBellard's formulaの両方で誤っていることから、共通の関数であるpow_modを調べました。
しばらくしてからk * k * nが64bit整数から溢れていることに気が付き、まず次のように書き直しました。
(((k * k) % d) * if m % 2 == 0 { 1 } else { n % d }) % d
つまりk * k * nを一気に計算せずに、一回掛けては剰余を取るようにしたのです。
これで10億桁目を計算してみました。
$ ./pihex_bbp 1000000000 1000000000: 73ccc138ae45 $ ./pihex_bellard 1000000000 1000000000: 5895585a7ef4
Bellard's formulaは5895585aまで合ってますが、BBP formulaの結果はまだ間違っています。
これではだめですね。
そこで であることを用いて次のように書き換えてみました。
let l = if k <= d / 2 { k } else { d - k }; (((l * l) % d) * if m % 2 == 0 { 1 } else { n }) % d
つまり と
の小さい方を二乗するというわけです。
実行結果は以下のようになります。
$ ./pihex_bbp 1000000000 1000000000: 58955859f966 $ ./pihex_bellard 1000000000 1000000000: 5895585a7ef4
ようやく両者のプログラム共に、10億桁目の正しい値を得られたようです。
pow_modの第三引数は8i + jですから、10億桁目を求める時は80億、つまり まで大きくなります。
上記のコードにおける
kはすなわちpow_modの返り値ですから、これも最大で くらいになります。
それを二乗すると
となり、unsigned 64bit整数の限界の
を簡単に超えてしまうのです。
を超えたら
を使うようにすると、掛け算の最大値を4分の1におさえることができます。
そうすると最大でも
つまりunsigned 64bit整数の限界をぎりぎりで下回り、オーバーフローすることなく計算できるのです。
100億桁目への挑戦
ところが100億桁となると、またいともたやすく64bit整数の限界を超えてしまいます。 あらためて状況を整理すると次のようになります。
が
を超える数で、
に対して、
を求めたい
の最大値の二乗が64bit整数の最大値を超えているので、単純な二乗はできません。
これまで64bit整数の限界と戦ってきたわけですが、よく調べてみるとRustでは128bit整数を使うことができます (記事執筆時点では#![feature(i128_type)]が必要ですが、そのうち書かなくてもよくなるでしょう*6 )。
つまり、最初のコード中のu64をu128に書き換えるだけで、オーバーフローの苦労なしに10億桁、そして100億桁を突破することができます。
以下の結果が、実際に128bit整数にして実行してみた様子です。
$ ./pihex_bbp 1000000000 1000000000: 58955859f966 # ./pihex_bbp 1000000000 7120.63s user 20.24s system 375% cpu 31:47.46 total $ ./pihex_bbp 10000000000 10000000000: 21c73c6889eb # ./pihex_bbp 10000000000 139050.65s user 365.69s system 373% cpu 10:21:40.48 total $ ./pihex_bellard 1000000000 1000000000: 5895585a7ef8 # ./pihex_bellard 1000000000 2643.48s user 6.98s system 374% cpu 11:47.56 total $ ./pihex_bellard 10000000000 10000000000: 21c73c65d3f4 # ./pihex_bellard 10000000000 80295.34s user 182.43s system 384% cpu 5:48:59.45 total
やった！
100億桁目を求めることができました！
実装上はdが大きくてオーバーフローの可能性がある場合にのみ128bitで計算するようにしています。
64bitで計算できる範囲のものを全て128bitで計算してしまうと、計算時間が余計にかかってしまうからです。
計算がうまく行ってしまうと感覚が麻痺してしまいますが、100億桁目の計算に成功したことに大きな喜びを感じています。
100億桁ですよ、100億桁！
公式をほぼそのまま書き下したコードを実行したら、ノートPCの上で現実的な時間内にこれだけの計算ができるというのは、本当に驚きました。
さらに計算機リソースの並列度を上げたり、計算結果をキャッシュできるものはキャッシュするとか、fract()の呼び出しを減らすなど、改善の余地はあると思います。
私が五年前に計算したのは円周率1億桁でした。
Bellard's formulaを使えば1億桁目のcb84を計算するのに30秒もかかりません。
10進数と16進数で見ている桁数が異なることに注意は必要ですが、十何分も掛けて計算した1億桁の最後の数字を30秒もかけずに検証できるということを肌で感じることができました。
数式を変形したり、コードを書いて試行錯誤する中で、円周率計算世界記録の技術にまた一歩近づけたのではないかなと思います。
おわりに
世界新記録の結果を検証するのは大変なことです。 同じ結果を持っている人がその人しかいないからです。 全く異なる方法で計算すると、元の計算以上の時間がかかってしまいます。 そういう状況で、特定の桁を直接求めることができるBBP formulaはとても重要な存在なのです。 浮動小数点を書き換えていくだけなので、巨大整数をファイルに書く必要もありませんし、メモリーも一定量しか消費しません。
円周率の16進数表現は、Blowfish暗号でも使われます。 Blowfish暗号の実装のためにBBP formulaを用いることもできますが、Blowfishで必要な桁数は決まっているので、ソースコードにべた書きされることが多いようです。 逆にBBPアルゴリズムを実装したら、検証データとしてBlowfish暗号の実装コードを参照することもできます。
浮動小数点演算に依存している以上、誤差は無視できるものではありません。
53bitの誤差が16進数表現の1桁をひっくり返すこともありえます。
ただし、ひっくり返ったかどうかは計算結果を見れば明らかに分かります。
その次の桁を見て000とかfffとなっているかどうかを見ればよいのです。
このアルゴリズムには、任意の桁の数字を正確に求めるよりも、小数点以下深い特定の桁を高速に求めることが要求されます。
間違える桁があることが否定できなくても、間違えたことがすぐに分かるのならそれでいいわけです。
本記事を書くにあたって、Rustを初めて真面目に書いてみました (以前、チュートリアルを少し触ったくらいでした)。 かなり書き心地がよくてびっくりしました。 普段書いている言語に対する不満が、いくつも解決されているように思えます。 今後はRustでもう少し大きなものを書いてみたいと思っています。
円周率計算は奥深い世界です。 円周率を計算する公式、Chudnovsky algorithmやBBP formula、Bellard's formulaなど、数論や解析学の結果によって何百日もの効率化につながることがあります。 円周率の16進数表現というわりとニッチな分野を調べていたのですが、Blowfish暗号と円周率の関係性を見つけるなど、意外な出会いもありました。
y-cruncherは円周率計算記録の多くで使われていることから分かるように、速度と安定性に評判があるソフトウェアです。 一方でy-cruncherのソースコードは公開されておらず、トータルで30万行のC++やインラインアセンブラで書かれているなど、なかなか謎の多いソフトウェアです。 ここまでのソフトウェアを書き上げた作者の情熱がどういうところにあるのか、とても気になります。
近年の円周率計算の世界記録は、y-cruncherの一強となっています。 これは面白くありませんので、いつかは兆桁台にも挑戦してみたいなと思っています。 みなさんも、ぜひ世界記録に挑戦してみてください。 Fabrice Bellard氏 *7 のComputation of 2700 billion decimal digits of Pi using a Desktop Computer (2009) は、円周率計算に関する素晴らしい資料です。 円周率計算の世界に少しでも興味を持っていただけたら幸いです。
*1:Chronology of computation of π https://en.wikipedia.org/wiki/Chronology_of_computation_of_%CF%80
*2:David H. Bailey, Peter Borwein, and Simon Plouffe, On the Rapid Computation of Various Polylogarithmic Constants, Mathematics of Computation 66, 903-913, 1997.
*3:Fabrice Bellard, A new formula to compute the n'th binary digit of pi, 1997.
*4:David V. Chudnovsky, Gregory V. Chudnovsky, The Computation of Classical Constants, Proceedings of the National Academy of Sciences of the United States of America 86 (21) 8178–8182, 1989.
*5:David H. Bailey. The BBP Algorithm for Pi, 2006.
*6:Tracking issue for 128-bit integer support (RFC 1504) https://github.com/rust-lang/rust/issues/35118
*7:Bellard's formulaの発見者であり、2009年12月に円周率を2.7兆桁計算して世界記録をとったプログラマです。FFmpegの作者でもあります。
