はじめに

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

\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 (The GNU Multiple Precision Arithmetic Library) 等、それ用のものを使ったほうが明らかに良いでしょう。

その他

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

.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)
任意精度演算についての紹介記事(Wikipedia)
よりちゃんとした多倍長演算を行いたいときは GMP
著名な円周率計算プログラム