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の実行結果をテキストエディタに貼り付けて、整形して、比較ツールで
比較してみました。
上図の左が、インターネットで探した1000桁の答えです。
右側が、上記のJavaで計算したπの値です。
大体500桁合っていますね。
バグってなくて、よかった、良かった。
おしまい。