円周率の小数点1000桁を出力したい

:::note info
以下の記事を自分なりに解釈し、自学した内容になります！
https://xn--w6q13e505b.jp/program/spigot.html
http://www.kk62526.server-shared.com/pi/Spigot.html
:::

# πをプログラムで計算しようと思い立った

3月14日なので計算しようと思い立ちました。
ただ、色々な近似式を見ても、項が増えていくとプログラムで計算できなくなる、、のをどうやって計算しているのだろうかというのは兼ねてからの疑問でした。

有名らしいライプニッツの級数は以下のような形で、arctanのマクローリン展開しており、これに1を代入することで${\pi}$の近似値が得られます。
(実はライプニッツより300年前にインドで発見されていたらしい)

```math
\arctan x = x-\frac{x^3} {3}+\frac{x^5} {5}-\cdots
```

```math
\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倍するを繰り返していく・・・みたいな感じです。

ネイピア数の例を用いて説明します。

```math
e = 1 + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + \frac{1}{4!} + \frac{1}{5!} + \frac{1}{6!} \cdots
```
これを式変形し、
```math
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桁目を計算します
```math
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]
```
さらに式変形し、
```math
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]
```
これを計算することで、
```math
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桁計算でき、出力できます。
```math
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]
```
これをもう一回やると
```math
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変換を行うことで
```math
\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
```
以下のように変換されます
```math
\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})
```
これを式変形すると

```math
\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$を代入し、式変形すると

```math
\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桁出力する。(数字の羅列でフォーマットはできてないです)

```rs
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桁でフォーマットするようにしたコードが以下

```rs
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)

```math
\frac{1}{\pi}=\frac{2\sqrt{2}}{99^2}\sum_{n=0}^\infty \frac{(4n)!(1103+26390n)}{(4^n99^nn!)^4}
```

Chudnovskyの級数(1988)

```math
\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のアルゴリズムは結構感動しました。無限に続くような数字の近似値ってどう出すんだろうというのはずっと長い間疑問には思っていたものの調べたこともなく、今回ついでに自分で手を動かして実装できたのは大変勉強になりました。