はじめに
この記事では、以下の式を使い、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しましたが
）
/**
* 階乗を計算します。
*/
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;
}