πをJavaで求めてみる（約５００桁）


## PI の求め方

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

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

## arctanの求め方

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

https://math-fun.net/20211121/20366/

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

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

https://ja.wikipedia.org/wiki/%E3%83%9E%E3%83%81%E3%83%B3%E3%81%AE%E5%85%AC%E5%BC%8F

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

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

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

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

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

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

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

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

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

```Java
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回繰り返して、終了条件とします。


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

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

http://www.koubou-kazu.co.jp/kzu/archi/sitanni/ensyuritu.html


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

![pi-kotae-awase.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1339875/5f2818ff-4953-d9bd-1d3d-008bd96b9500.png)


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

大体５００桁合っていますね。

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