以下の記事を自分なりに解釈し、自学した内容になります！
https://xn--w6q13e505b.jp/program/spigot.html
http://www.kk62526.server-shared.com/pi/Spigot.html
πをプログラムで計算しようと思い立った
3月14日なので計算しようと思い立ちました。
ただ、色々な近似式を見ても、項が増えていくとプログラムで計算できなくなる、、のをどうやって計算しているのだろうかというのは兼ねてからの疑問でした。
有名らしいライプニッツの級数は以下のような形で、arctanのマクローリン展開しており、これに1を代入することで${\pi}$の近似値が得られます。
(実はライプニッツより300年前にインドで発見されていたらしい)
\arctan x = x-\frac{x^3} {3}+\frac{x^5} {5}-\cdots
\arctan{1}=\frac{\pi}{4}=\frac{1}{1}-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\frac{1}{11}+\frac{1}{13}\cdots
といってもこれをそのまま実装するとプログラムだと割とすぐに限界がきそうですよね。。
Spigotアルゴリズム
このアイデアで、無限に精度がなくても計算できるようになります。
詳しくはこの後書きますがざっくりいうと1桁出力し、管理している分子を10倍するを繰り返していく・・・みたいな感じです。
ネイピア数の例を用いて説明します。
e = 1 + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + \frac{1}{4!} + \frac{1}{5!} + \frac{1}{6!} \cdots
これを式変形し、
e = 2 + \frac{1}{2} \bigg( 1 + \frac{1}{3} \bigg( 1 + \frac{1}{4} \bigg( 1 + \frac{1}{5} \bigg( 1 + \frac{1}{6} \bigg)\bigg)\bigg)\bigg)\bigg)
2桁目を計算します
e = 2 + \frac{1}{10} \bigg[ 10 \times \frac{1}{2} \bigg( 1 + \frac{1}{3} \bigg( 1 + \frac{1}{4} \bigg( 1 + \frac{1}{5} \bigg( 1 + \frac{1}{6} \bigg)\bigg)\bigg)\bigg)\bigg) \bigg]
さらに式変形し、
e = 2 + \frac{1}{10} \bigg[\frac{1}{2} \bigg( 10 + \frac{1}{3} \bigg( 10 + \frac{1}{4} \bigg( 10 + \frac{1}{5} \bigg( 10 + \frac{10}{6} \bigg)\bigg)\bigg)\bigg)\bigg) \bigg]
これを計算することで、
e = 2 + \frac{1}{10} \bigg[7 + \frac{1}{2} \bigg( 0 + \frac{1}{3} \bigg( 1 + \frac{1}{4} \bigg( 0 + \frac{1}{5} \bigg( 1 + \frac{4}{6} \bigg)\bigg)\bigg)\bigg)\bigg) \bigg]
1桁計算でき、出力できます。
e = 2.7 + \frac{1}{10} \bigg[\frac{1}{2} \bigg( 0 + \frac{1}{3} \bigg( 1 + \frac{1}{4} \bigg( 0 + \frac{1}{5} \bigg( 1 + \frac{4}{6} \bigg)\bigg)\bigg)\bigg)\bigg) \bigg]
これをもう一回やると
e = 2.71 + \frac{1}{100} \bigg[\frac{1}{2} \bigg( 1 + \frac{1}{3} \bigg( 1 + \frac{1}{4} \bigg( 3 + \frac{1}{5} \bigg( 1 + \frac{4}{6} \bigg)\bigg)\bigg)\bigg)\bigg) \bigg]
と計算できます。これなら分数の計算もいらずに無限の精度がなくても計算できる・・・！！！すごい！
Eulerのarctan展開
Spigotのアルゴリズムを適用すべく、arctanの級数を変換します。
以下のEuler変換を行うことで
\sum_{n=0}^{\infty}(-1)^n a_n x^n = \frac{1}{1+x} \sum_{n=0}^{\infty} (-1)^n \left(\frac{x}{1+x}\right)^n \Delta^na_0
以下のように変換されます
\arctan x = \frac{y}{x} \bigg(1+\dfrac{2}{3}y + \dfrac{2\cdot4}{3\cdot5}
y^2 + \dfrac{2\cdot4\cdot6}{3\cdot5\cdot7}y^3 + \cdots\bigg)
\quad ({\rm ただし\ }y = \dfrac{x^2}{1+x^2})
これを式変形すると
\arctan x = \dfrac{y}{x} \left(1 + \dfrac{2}{3}y
\left(1 + \dfrac{4}{5}y \left(1 + \dfrac{6}{7}y
\left(1 + \cdots \right) \right) \right) \cdots \right)
$x=1$を代入し、式変形すると
\pi = 2+\frac13\left(2+\frac25\left(2+\frac37\left(2+\cdots\left(2+\frac{k}{2k+1}\left(2+\cdots\right) \cdots \right.\right.\right.\right)
が得られる。これを4桁ずつ出力するプログラムを作成する。
実コードたち(Rustです)
小数点以下2399桁出力する。(数字の羅列でフォーマットはできてないです)
fn main() {
let base = 10000;
let mut n = 8400;
let mut numerator = vec![2000; n];
let mut out = 0;
while 0 < n {
let mut temp = 0;
for i in (1..n).rev() {
let denom = 2 * i - 1;
temp = temp * i + numerator[i] * base;
numerator[i] = temp % denom;
temp /= denom;
}
print!("{:<04}", out + temp / base);
out = temp % base;
n -= 14;
}
}
小数点以下1000桁でフォーマットするようにしたコードが以下
fn main() {
let base = 10000;
let mut n = 8400;
let mut numerator = vec![2000; n];
let mut out = 0;
let mut ans = String::new();
for _ in 0..251 {
let mut temp = 0;
for i in (1..n).rev() {
let denom = 2 * i - 1;
temp = temp * i + numerator[i] * base;
numerator[i] = temp % denom;
temp /= denom;
}
let print = format!("{:<04}", out + temp / base);
ans.push_str(&print);
out = temp % base;
n -= 14;
}
ans.remove(0);
println!("3.{}", ans.chars().take(1000).collect::<String>());
}
code golf用に314byteに圧縮
fn main(){let b=10000;let mut n=8400;let mut m=vec![2000;n];let mut o=0;let mut a=String::new();for _ in 0..251{let mut t=0;for i in (1..n).rev(){let d=2*i-1;t=t*i+m[i]*b;m[i]=t%d;t/=d;}let s=format!("{:<04}",o+t/b);a.push_str(&s);o=t%b;n-=14;}a.remove(0);println!("3.{}",a.chars().take(1000).collect::<String>())}
その後なんやかんやして291byteに改善
fn main(){let b=10000;let mut n=8400;let mut v=vec![b/5;n];let mut a=String::new();let mut c=0;while 0<n{c%=b;let d=c;for i in(1..n).rev(){let e=2*i-1;c=c*i+b*v[i];v[i]=c%e;c/=e;}a.push_str(&format!("{:<04}",d+c/b));n-=14;}a.remove(0);print!("3.{}",a.chars().take(1000).collect::<String>())}
おまけ
自分が実装した関数は$O(n^2)$で動作しますが、現代はもっとはやく計算できるようになっているらしいです。
収束が速い有名な級数たちを紹介します。Binary Splittingというテクニックを使うと$O(n(\log n)^3)$で計算できるらしいです。項数が$2^n$の場合に分割統治することで計算量が落ちるらしいです。(そこまで理解できていないです)
現代のすごい桁数を計算しているものはChudnovskyの級数(Ramanujanの級数を改良したもの)を用いて計算しており、1つの項で14桁ほど求まるらしいです。現在は100兆桁以上円周率は求まっているのだとか。
Ramanujanの級数(1914)
\frac{1}{\pi}=\frac{2\sqrt{2}}{99^2}\sum_{n=0}^\infty \frac{(4n)!(1103+26390n)}{(4^n99^nn!)^4}
Chudnovskyの級数(1988)
\frac{1}{\pi}=12\sum_{n=0}^\infty \frac{{(-1)}^{n}(6n)!}{(3n)!{(n!)}^3} \frac{13591409+545140134n}{{640320}^{3n+\frac{3}{2}}}
昨今の桁数世界記録のアルゴリズム自体は2009年に発表されたy-cruncherの独壇場となっているみたいです。並列化とかマシンパワーで桁数が伸びているみたいです。
感想
数学っぽい知識は理解するの時間かかってしまいまして、大人になってから思うともうちょっと勉強していればよかったなぁなんて思いました。
それとは別に、Spigotのアルゴリズムは結構感動しました。無限に続くような数字の近似値ってどう出すんだろうというのはずっと長い間疑問には思っていたものの調べたこともなく、今回ついでに自分で手を動かして実装できたのは大変勉強になりました。