Javaで円周率

# はじめに

この記事では、以下の式を使い、Javaで円周率を計算してみます。

```math
\frac{4}{\pi}=\sum_{n=0}^\infty \frac{(-1)^n(4n)!(1123+21460n)}{882^{2n+1}(4^nn!)^4}
```

ラマヌジャン(1887/12/22 - 1920/4/26)の公式です。
しげしげ見ていて、どうやってこれを思いついたのかとか、ほんとうにこれで円周率が出てくるのかとか、不思議な気分になります。

（どの数式に関するやり取りか、自分は知らないのですが）以下のやり取りもミステリアスです。

ハーディ教授「どうやってこの数式を思いついたのか教えて欲しい」
ラマヌジャン「どうか信じてほしい、ナマギーリ女神が枕元に立って教えてくれたんだ」

実際に計算してみたくなります。

# プログラム

基本的には、工夫のないプログラムです。

- 分子と分母はそれぞれBigIntegerで計算しています。
- 分数を割り算するところからBigDecimalで計算します(このプログラムでは、ここではじめて精度指定します)。
- 階乗は同じ計算を何度もするのはあまりにひどいので、前の計算結果を再利用するよう「工夫」しています。
- 最後の方の桁は誤差を含みます（誤差は「繰り返し回数」に比例します）。

``` PamaPi.java
/**
 * 円周率を計算します。
 */

import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.RoundingMode;

public class RamaPi {
	/**
	 * 階乗を計算します(直前の計算結果を再利用します)。
	 */
	class Facto {
		/** 直前の引数 */
		int n = 1;
		/** 直前の値 */
		BigInteger val = BigInteger.ONE;

		BigInteger calc(final int arg) {
			for ( ; n <= arg ; n++) {
				val = val.multiply(BigInteger.valueOf(n));
			}

			return val;
		}
	}

	/** 階乗計算用１ */
	Facto workFacto1 = new Facto();
	/** 階乗計算用２ */
	Facto workFacto2 = new Facto();
	/** 定数「882」 */
	static final BigInteger VAL_882 = BigInteger.valueOf(882);
	/** 定数「4」 */
	static final BigInteger VAL_4 = BigInteger.valueOf(4);

	/**
	 * ラマヌジャンの公式の分子を計算します。
	 */
	BigInteger bunsi(final int n) {
		final int sign = (n&1) == 0 ? 1 : -1;
		final BigInteger a = workFacto1.calc(4 * n);
		final BigInteger b = BigInteger.valueOf(sign * (1123 + 21460*n));

		return a.multiply(b);
	}

	/**
	 * ラマヌジャンの公式の分母を計算します。
	 */
	BigInteger bunbo(final int n) {
		final BigInteger a = VAL_882.pow(2*n+1);
		final BigInteger b = VAL_4.pow(n);
		final BigInteger c = workFacto2.calc(n);
		final BigInteger d = b.multiply(c).pow(4);
		return a.multiply(d);
	}

	/**
	 * ラマヌジャンの公式から、円周率を求めます。
	 */
	void calc(final int digits) {
		BigDecimal sum = BigDecimal.ZERO;

		/** 加算を打ち切る桁数（この桁数未満は足せません） */
		final BigDecimal thr = BigDecimal.ONE.divide(BigDecimal.TEN.pow(digits));

		/** 計算開始した時間 */
		long t0 = System.nanoTime();

		int iter = 0;
		for ( ;; iter++ ) {
			final BigDecimal bunbo = new BigDecimal(bunbo(iter));
			final BigDecimal bunsi = new BigDecimal(bunsi(iter));

			final BigDecimal adder = bunsi.divide(bunbo, digits, RoundingMode.HALF_EVEN);

			if (adder.abs().compareTo(thr) < 0) {
				break;
			}

			/* 加算の結果で、sumを書き換えます */
			sum = sum.add(adder);

			/* 桁数が小さい場合、収束の様子を表示します */
			if (digits < 100) {
				System.out.println( BigDecimal.valueOf(4).divide(sum, digits, RoundingMode.HALF_EVEN));
			}
		}

		/** 計算終了した時間 */
		long t1 = System.nanoTime();
		System.err.println("繰り返した回数: " + iter);
		double elapsed = (t1 - t0) / ((double)1000 * 1000 * 1000);
		System.err.println("所要した秒数: " + String.format("%.6f", elapsed));

		System.out.println(BigDecimal.valueOf(4).divide(sum, digits, RoundingMode.HALF_EVEN));
	}

	static public void main(String arg[]) {
		int digits = Integer.parseInt(arg[0]);

		new RamaPi().calc(digits);
	}
}

```

# 実行例

引数で、桁数を指定します。
最後の数桁に誤差が乗る可能性があります（誤差は「繰り返した回数」に比例します）。

```
java RamaPi 50
3.14158504007123775601068566340160284951024042742654
3.14159265359762201426827517920229156338712376594028
3.14159265358979322988767708881201158250039662639857
3.14159265358979323846265313702144396820307036308333
3.14159265358979323846264338326814215966511742779251
3.14159265358979323846264338327950289764449222972749
3.14159265358979323846264338327950288419715329619277
3.14159265358979323846264338327950288419716939939455
3.14159265358979323846264338327950288419716939937511
繰り返した回数: 9
所要した秒数: 0.003002
3.14159265358979323846264338327950288419716939937511
```

# 所要時間について

手元の環境で、少し測りました。
なかなか時間がかかりますが、Javaの多倍長計算を速いと捉えるか、遅いと捉えるか。

|    桁数    |    繰り返した回数    |    所要した秒数    |    1回あたりμs    |
|:-----------|---------------------:|:------------------:|:------------------:|
|           10    |    2    |    0.001129            |      564.5 |
|           100    |    18    |    0.003077          |      170.9 |
|           1000    |    170    |    0.049033        |      288.4 |
|           10000    |    1698    |    4.237159      |     2495.3 |
|           20000    |    3395    |    19.83075      |     5841.1 |
|           40000    |    6790    |    118.651706    |    17474.4 |

なお、同じ環境で「Superπ」を動かしてみると、「838万桁」が「152秒」でした（桁数で200倍、違います）。

より速度を追求したい場合、多倍長演算としては、Java の BigInteger / BigDecimal ではなく [GMP](https://gmplib.org/) (The GNU Multiple Precision Arithmetic Library) 等、それ用のものを使ったほうが明らかに良いでしょう。

# その他

当初、素朴な階乗の実装をしていましたが、たったの 2 万桁で 113 秒も所要しました。
（一番最初はとりあえず再帰で書いてみたところ、すぐstack overflowしましたが:hushed:）

``` .java
	/**
	 * 階乗を計算します。
	 */
	static final BigInteger facto(final int a) {
		BigInteger work = BigInteger.ONE;

		for (int i=2 ; i<=a ; i++) {
			work = work.multiply(BigInteger.valueOf(i));
		}
		return work;
	}
```

# 参考

* [シュリニヴァーサ・ラマヌジャン(Wikipedia)](https://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A5%E3%83%AA%E3%83%8B%E3%83%B4%E3%82%A1%E3%83%BC%E3%82%B5%E3%83%BB%E3%83%A9%E3%83%9E%E3%83%8C%E3%82%B8%E3%83%A3%E3%83%B3)

* [任意精度演算](https://ja.wikipedia.org/wiki/%E4%BB%BB%E6%84%8F%E7%B2%BE%E5%BA%A6%E6%BC%94%E7%AE%97) についての紹介記事(Wikipedia)

* よりちゃんとした多倍長演算を行いたいときは  [GMP](https://gmplib.org/)

* [著名な円周率計算プログラム](http://円周率.jp/program/introduction.html)