ゆっちのBlog

円周率を Machin’s formula で求める

Java

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

今回はマチンの公式 ( 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;
    }
}

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<BigDecimal> {

    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<Future<BigDecimal>> 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;
    }
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

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<BigDecimal> {

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<Future<BigDecimal>> 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 桁で確認したところ問題なさそうでした。

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

Hatena タグ: Java

TAGS: Java | 2013年5月6日7:23 PM

« もっと Lambda その７もっと Lambda その８ »

Comment

lic aao exam syllabussays:
2016年3月2日 5:46 PM

Good post. I’m going through some of these issues as well..

返信
cream pembesar payudarasays:
2017年6月17日 7:11 PM

Hello to every body, it’s my first pay a visit of this webpage; this web site carries amazing
and actually excellent information in favor of readers.

返信

円周率を Machin’s formula で求める

Comment

Leave a Reply (or Cancel)

Trackback URL

Recent Comments

Java

NetBeans

カテゴリー

アーカイブ

RSS FEED

ねこび～ん

メタ情報

ゆっちのカレンダー

2024年11月
月	火	水	木	金	土	日
« 1月
				1	2	3
4	5	6	7	8	9	10
11	12	13	14	15	16	17
18	19	20	21	22	23	24
25	26	27	28	29	30

ゆっちのBlog

円周率を Machin’s formula で求める

Comment

Leave a Reply (or Cancel)

Trackback URL

Popular post ranking

Recent Comments

Java

NetBeans

カテゴリー

アーカイブ

RSS FEED

ねこび～ん

メタ情報

ゆっちのカレンダー