2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

πをJavaで求めてみる(約500桁)

Posted at

PI の求め方

tan(45度) = tan(π / 4) =1 は、ご存知かと思います。

tanの逆関数を arctanとすると
π / 4 = arctan(1) となります。
両辺を4倍すると
π = 4 * arctan(1) となります。

arctanの求め方

テーラー展開により求める事ができます。
私は、こちらのページを参照しました。

ところが、arctan(1)は収束が遅い

WikiPediaによると、色々な公式が考案されています。

マチンの公式は有名ですが、これでもまだまだ収束が遅いので、このページの下の方にある公式を、あさります。

3項以上の公式のうち、一番左の計算が遅いので、
これが小さな値のものを探します。

ガウスによる公式 ==> 一番左が (1 / 18)

ストーマーの公式(1) ==> 一番左が (1 / 57)

ストーマーの公式(2) ==> 一番左が (1 / 8)

高野喜久雄による公式 ==> 一番左が (1 / 49)

シムソンの公式 ==> 一番左が (1 / 10)

となっていますので、
ストーマーの公式(1)が最も収束が早く、
高野喜久雄の公式が2番目に収束が早そう、と予想できます。

素直に ストーマーの公式(1)を実装

package calcpi;

import java.math.BigDecimal;
import java.math.MathContext;

public class CalcPI {

    static MathContext mi = new MathContext(512);  // (1)

    static BigDecimal limit; // (2)
    /*
    *   arctan の求め方
    *   https://math-fun.net/20211121/20366/#google_vignette
    */
    static BigDecimal arctan(BigDecimal arg) {

        BigDecimal accum = arg;
        BigDecimal bunshi = arg;
        BigDecimal square = arg.multiply(arg, mi);
        BigDecimal bunbo = new BigDecimal("3");
        BigDecimal two = new BigDecimal("2");

        boolean minus = true;

        for (int i = 0;; i++) {

            bunshi = bunshi.multiply(square, mi);
            
            BigDecimal kou = bunshi.divide(bunbo, mi);
            
            System.out.println("" + i +" = " + kou);

            if (kou.compareTo(limit) < 0) {  // (3)
                break;                       // (4)
            }

            if (minus) {
                accum = accum.subtract(kou, mi);
            } else {
                accum = accum.add(kou, mi);
            }
            bunbo = bunbo.add(two, mi);
            minus = !minus;
        }
        return accum;
    }

    /*
     * 
     * https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%81%E3%83%B3%E3%81%AE%E5%85%AC%E5%BC%8F
     * マチンの公式
     */
    public static void main(String[] args) {

        BigDecimal stop = new BigDecimal("1");
        
        for (int i = 0; i < 512; i++) {
            stop = stop.divide(new BigDecimal("10"), mi);
        }
        limit = stop;   // (5)
        
        
        // ストーマー (Fredrik Carl Mulertz Stormer, 1874-1957) による公式 
        
        BigDecimal one = new BigDecimal("1");

        BigDecimal d57 = one.divide(new BigDecimal("57"), mi);
        BigDecimal d239 = one.divide(new BigDecimal("239"), mi);
        BigDecimal d682 = one.divide(new BigDecimal("682"), mi);
        BigDecimal d12943 = one.divide(new BigDecimal("12943"), mi);
        
        BigDecimal atan57 = arctan(d57);
        BigDecimal atan239 = arctan(d239);
        BigDecimal atan682 = arctan(d682);
        BigDecimal atan12943 = arctan(d12943);



        BigDecimal kou1 = atan57.multiply(new BigDecimal("44"), mi);
        BigDecimal kou2 = atan239.multiply(new BigDecimal("7"), mi);
        BigDecimal kou3 = atan682.multiply(new BigDecimal("12"), mi);
        BigDecimal kou4 = atan12943.multiply(new BigDecimal("24"), mi);
        
        BigDecimal pi4 = kou1.add(kou2, mi).subtract(kou3, mi).add(kou4, mi);
        
        BigDecimal ans = pi4.multiply(new BigDecimal("4"), mi);
        System.out.println("ans=" + ans);
    }

}


(1)は、だいたい512桁ぐらい求めてね、という指示です。
 (細かい事は、私もよくわかっていません)

(2)は、500桁ほど求まったら、arctanの計算をやめるための数です。

(3), (4)のところで、テーラー展開の項が、(2)の数値以下になったら、ループを脱出します。
これが、arctanの答えになります。

mainの最初から、(5)のあたりが、(2)の変数を初期化しています。
1➗ 10を512回繰り返して、終了条件とします。

答え合わせのページを探す

以下のページに1000桁のっていました。

Javaの実行結果をテキストエディタに貼り付けて、整形して、比較ツールで
比較してみました。

pi-kotae-awase.png

上図の左が、インターネットで探した1000桁の答えです。
右側が、上記のJavaで計算したπの値です。

大体500桁合っていますね。

バグってなくて、よかった、良かった。
おしまい。

2
0
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?