一年近く前に円周率をモンテカルロ法を利用して求めるプログラムを組んだ。

今回は マチンの公式 ( Machin’s formula ) を利用してみようと思う。

そして桁数と精度も上げてみたいと思ってこの方法を選択してみた。

まず、マチンの公式とは次のようになっている。

この公式の証明は三角関数の公式を利用して証明できます。

倍角公式を二度使い、加法定理でとどめを刺して終わりです。

詳しくはググってみてください。

それではプログラムをこれにそって組んでみましょう。

プログラムを組むにあたって都合の良いように式を変えてあります。( 6 行目 )

このプログラムでは arcTan(double x) の n 番目 (i + i + 1) の項が 0 になるまで計算をします。

十分小さな値をとるまでに変更する場合は、if 文の条件は絶対値をとるように変更すればいいでしょう。

package jp.yucchi.highprecisionpi;

public class HighPrecisionPi {

    public static void main(String[] args) {
        double pi = 16.0 * arcTan(1.0 / 5.0) - 4.0 * arcTan(1.0 / 239.0);
        System.out.println(pi);
    }

    private static double arcTan(double x) {
        double ret = 0.0;
        double sig = x;
        double sqx = x * x;
        for (int i = 0; sig != 0.0; i++) {
            ret += sig / (double) (i + i + 1);
            sig = -sig * sqx;
        }
        return ret;
    }
}

さて、このプログラムの実行結果を確認してみましょう。

3.141592653589794

最後の桁が 3 であるところが 4 になってます。

モンテカルロ法より速くて精度は高いですね。

でも、もっと桁数をとりたくないでしょうか？

そこで、このプログラムを次のように変更してみました。

百万桁の計算をさせるために BigDecimal を使ってみました。

基本的に先ほどのプログラムと同じです。

ちょっと処理時間が気になるので時間計測をいれてます。

と言っても百万桁も計算させると凄く時間がかかってしまうのは解りきっていることなので一万桁で試してみます。

普段使う機会が少ない BigDecimal を使ってますのでちょっとおかしいところがあるかもしれませんが気にしないでね。(^_^;)

あと、桁数を多くとるので 0 収束を使うと非常に時間がかかってしまうので誤差がでないように適当にループさせてます。

それではいつものように未熟で汚いコードですがご覧ください。

package jp.yucchi.highprecisionpi;

import java.lang.management.ManagementFactory;
import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.*;
import java.util.logging.Level;
import java.util.logging.Logger;

public class HighPrecisionPi implements Callable {

    private static final int KETA = 10_000;
    private static final int LIMIT_NUMBER = KETA + 2;
    private static final int LIMIT_NUMBER2 = KETA ;
    private static long startTime;
    private static long time;
    private final BigDecimal s_q;
    private final BigDecimal _limit;
    private BigDecimal ss;

    public HighPrecisionPi(BigDecimal _limit, BigDecimal ss, BigDecimal s_q) {
        this._limit = _limit;
        this.ss = ss;
        this.s_q = s_q;
    }

    public static void main(String[] args) {
        System.out.println("<--- 円周率計算開始 --->");
        startTime = System.nanoTime();
        //精度KETA桁 四捨五入で丸める
        int setPrecision = KETA;
        RoundingMode setRoundingMode = RoundingMode.HALF_UP;
        MathContext mc2 = new MathContext(1, setRoundingMode);
        MathContext mc = new MathContext(setPrecision, setRoundingMode);
        BigDecimal f = BigDecimal.ZERO;  // 0
        BigDecimal r = BigDecimal.ZERO;  // 0
        BigDecimal one = BigDecimal.ONE;  // 1     
        BigDecimal two = new BigDecimal(2, mc); //2      
        BigDecimal n239 = new BigDecimal(239, mc); // 239
        BigDecimal sq239 = new BigDecimal(57121, mc); // 239^2

        int procs = ManagementFactory.getOperatingSystemMXBean().getAvailableProcessors();
        if (procs > 1) {
            procs = 2;
        }
        ExecutorService executor = Executors.newFixedThreadPool(procs);

        List> futures = new ArrayList<>();
        for (int i = 0; i < 2; i++) {
            final BigDecimal _limit;
            final BigDecimal ss;
            final BigDecimal s_q;
            if (i > 0) {
                _limit = new BigDecimal(LIMIT_NUMBER, mc);
                ss = new BigDecimal(0.2, mc2); // 1 / 5
                s_q = new BigDecimal(0.04, mc2); // 1 / 25
            } else {
                _limit = new BigDecimal(LIMIT_NUMBER2, mc);
                ss = one.divide(n239, mc); // 1 / 239        
                s_q = one.divide(sq239, mc); // 1 / 57121     
            }
            futures.add(executor.submit(new HighPrecisionPi(_limit, ss, s_q)));
        }

        try {
            f = futures.get(1).get().multiply(new BigDecimal(16, mc));
        } catch (InterruptedException | ExecutionException ex) {
            Logger.getLogger(HighPrecisionPi.class.getName()).log(Level.SEVERE, null, ex);
        }

        try {
            r = futures.get(0).get().multiply(new BigDecimal(4, mc));
        } catch (InterruptedException | ExecutionException ex) {
            Logger.getLogger(HighPrecisionPi.class.getName()).log(Level.SEVERE, null, ex);
        }

        System.out.println("\n" + KETA + "桁の円周率は、\n" + f.subtract(r, mc));
        time = System.nanoTime() - startTime;
        System.out.println("\n" + KETA + "桁の円周率計算時間は、"
                + (int) (time * 1e-9) / 3_600 + "時間"
                + (int) ((time * 1e-9) / 60) % 60 + "分"
                + (int) (time * 1e-9 % 60) + "秒"
                + Double.toString((time * 1e-9 % 60) % 1).substring(2) + "\n");

        System.out.println("<--- 円周率計算終了 --->");

        executor.shutdown();
    }

    @Override
    public BigDecimal call() throws Exception {
        int setPrecision = KETA;
        RoundingMode setRoundingMode = RoundingMode.HALF_UP;
        MathContext mc = new MathContext(setPrecision, setRoundingMode);
        MathContext mc2 = new MathContext(1, setRoundingMode);
        BigDecimal fr = BigDecimal.ZERO;  // 0
        BigDecimal p2 = BigDecimal.ZERO;  // 0
        BigDecimal one = BigDecimal.ONE;  // 1     
        BigDecimal two = new BigDecimal(2, mc); //2 

        for (int i = 0; p2.compareTo(_limit) < 0; i++) {
            BigDecimal bi = new BigDecimal(i, mc);
            fr = fr.add(ss.divide((bi.multiply(two).add(one)), mc));
            ss = ss.multiply(s_q, mc).negate();
            p2 = p2.add(one);
        }
        return fr;
    }
}

少しでも計算時間が短くなればいいなぁって思って悪あがきしてます。(￣。￣;)

基本的な考え方は double 仕様のものと一緒で悪あがきの並行処理を取り入れただけです。

実行結果の抜粋を載せておきます。

<--- 円周率計算開始 --->

10000桁の円周率は、
3.1415926535897932384626433832795028841971693993751058209749445923078164062

略

4988727584610126483699989225695968815920560010165525637568

10000桁の円周率計算時間は、0時間3分22秒4562651430000244

<--- 円周率計算終了 --->

ちゃんと計算できてますね。

一番最後は丸めによって四捨五入されてます。

もっと効率よくプログラムを組むには収束値を設定したほうがいいのかなぁ・・・？

でも、速くても計算が甘くなるのはいやだ！

とりあえず昨年作ったモンテカルロ法で円周率を求めるプログラムより遙かに良いのは確かだ！

あと気のせいかもしれないが BigDecimal の除算って遅く感じる。

ゴールデンウィークの最後の休日にこんなことしているってのは平和でいいね。(^_^)

平成25年5月8日(水)　追記　計算のループの回数を決めている LIMIT_NUMBER の値を変更しました。

10000 桁で確認したところ問題なさそうでした。

ただし、それより大きな桁は確認していませんのであしからず。(；´Д｀)