Help us understand the problem. What is going on with this article?

Javaで円周率

More than 1 year has passed since last update.

はじめに

この記事では、以下の式を使い、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;
        }
    }

    /** 階乗計算用1 */
    Facto workFacto1 = new Facto();
    /** 階乗計算用2 */
    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しましたが:hushed:

    /**
     * 階乗を計算します。
     */
    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;
    }

参考

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away