円周率の計算式とその収束速度

はじめに

円周率を計算する式は数多く提案されている。どの式が、どの程度速いのか、つまり、収束の具合について、エクセルやprocessingで調べてみたい。

※現状で桁数の記録を狙うには、大量のメモリ、大量のディスク、計算効率、計算時間が必要で、資本がまぁまぁ必要となる。

思うこと

円周率πの定義は単純明快。だが無限に続く。そして、それを計算するための式がたくさんある。唯一の答えに対して、たくさんのアプローチが存在する。それが不思議。

BigDecimalの使用例(java)

import java.math.*;

MathContext mc = new MathContext(100);

BigDecimal s = new BigDecimal("2");
s = s.sqrt(mc);
println(s);

オイラー公式2乗の収束速度

\frac{1}{1^2}+
\frac{1}{2^2}+
\frac{1}{3^2}+
\frac{1}{4^2}+
\frac{1}{5^2}+
\cdots
=
\frac{\pi^2}{6}

__10項で誤差0.1
_100項で誤差0.01
1000項で誤差0.001
...

オイラー公式4乗の収束速度

\frac{1}{1^4}+
\frac{1}{2^4}+
\frac{1}{3^4}+
\frac{1}{4^4}+
\frac{1}{5^4}+
\cdots
=
\frac{\pi^4}{90}

_15項で誤差0.00001
_32項で誤差0.000001
_69項で誤差0.0000001
149項で誤差0.00000001
...

オイラー公式6乗の収束速度

\frac{1}{1^6}+
\frac{1}{2^6}+
\frac{1}{3^6}+
\frac{1}{4^6}+
\frac{1}{5^6}+
\cdots
=
\frac{\pi^6}{945}

__5項で誤差0.00001
__7項で誤差0.000001
_11項で誤差0.0000001
_18項で誤差0.00000001
_29項で誤差0.000000001
_46項で誤差0.0000000001
_72項で誤差0.00000000001
116項で誤差0.000000000001
...

オイラー公式8乗の収束速度

\frac{1}{1^8}+
\frac{1}{2^8}+
\frac{1}{3^8}+
\frac{1}{4^8}+
\frac{1}{5^8}+
\cdots
=
\frac{\pi^8}{9450}

__3項で誤差0.0001
__4項で誤差0.00001
__5項で誤差0.000001
__8項で誤差0.0000001
_10項で誤差0.00000001
_15項で誤差0.000000001
_20項で誤差0.0000000001
_28項で誤差0.00000000001
_39項で誤差0.000000000001
_55項で誤差0.0000000000001
_76項で誤差0.00000000000001
Doubleの限界

オイラー公式？3乗の収束速度

\frac{1}{1^3}-
\frac{1}{3^3}+
\frac{1}{5^3}-
\frac{1}{7^3}+
\frac{1}{9^3}-
\cdots
=
\frac{\pi^3}{32}

__4項で誤差0.001
__9項で誤差0.0001
_19項で誤差0.00001
_40項で誤差0.000001
_86項で誤差0.0000001
185項で誤差0.00000001

オイラー積（素数）の公式

\Big(\frac{1}{1-\frac{1}{2^2}}\Big)
\Big(\frac{1}{1-\frac{1}{3^2}}\Big)
\Big(\frac{1}{1-\frac{1}{5^2}}\Big)
\Big(\frac{1}{1-\frac{1}{7^2}}\Big)
\Big(\frac{1}{1-\frac{1}{11^2}}\Big)
\Big(\frac{1}{1-\frac{1}{13^2}}\Big)
\cdots
=
\frac{\pi^2}{6}

__3項で誤差0.1
__12項で誤差0.01
__54項で誤差0.001

マティン（マチン）

\pi=16\tan^{-1}\Big(\frac{1}{5}\Big)-4\tan^{-1}\Big(\frac{1}{239}\Big)

ガウス

\pi=48\tan^{-1}\Big(\frac{1}{18}\Big)+32\tan^{-1}\Big(\frac{1}{57}\Big)-20\tan^{-1}\Big(\frac{1}{239}\Big)

Stormer

\pi=
24\tan^{-1}\Big(\frac{1}{8}\Big)+
8\tan^{-1}\Big(\frac{1}{57}\Big)+
4\tan^{-1}\Big(\frac{1}{239}\Big)

arctanの近似式

以上の3式では、円周率の精度はarctanの精度に依存するため、arctanの近似式の収束速度について調べる。

arctan(x)について、近似した場合、０付近の精度は高いが、1もしくは、-1に近づくほど、急激に精度が悪くなることがグラフから読み取れる。よって、できるだけ小さい（０に近い）数を使うことが求められる。

例えば、マティンの式では、1/5や1/239といった、小さい数字が登場する。精度が高められやすい式に工夫されている。
以後、計算上、arctan(x)よりも、arctan(1/x)が登場する。

「マティン」の式で1000桁の円周率を計算するプログラムは以下の通り。

import java.math.*;

MathContext mc = new MathContext(1000);

BigDecimal _5 = new BigDecimal("0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441126281180736913601044564798867942393557475654952163032700522107470015645015560061286185526633257318692806643896806189528405825931124251613297313993397113233537821796084176648310525473039665725650488878155309384290579311695934192851806364919697519401708560949527368673738508400812367856158009329822514023246675549211026704574378815474839079978985020075223696837961392278354193255722328413846477441352909705465122438302697560518377574220877835853152464749330914587633823112490332030126805100670223312575050942448460267162254894079226140467995062365969282873058287872053603034570766066681374312566267431408992605741703545394046813622411076808100362137599259589071666485145255706799548810043132954668389229036883093278538312521768924841625558744394478718875141137978430670810332362273330875483507230420969155374569324447556530899212193854144645083686720518519002662277446939325408961915860225494493742927");
BigDecimal _239 = new BigDecimal("0.004184076002074723864538214959285452741048065307631950827019612887181778341422893273782605813622909454975450666444863756052458394789311865058922128833092800846271962330773375947634603318473414570331986015454814805992449830211460391253949527607796881558881273397853346518045742548135867464475197910232830977002064652827634653296910481838654356078919591451232220944636862766155208316796426465746551103251034352628244512693556704996844452479043317728393070863140193869519503705864107708558554045223553881423767708365156918252702002293089544950043585440934496440142418724950922838623954553335651171973747020234947597790974695011188854766739795731537093032782113089842583083677190910083909851655104192241678092053268649162667402716844424477315796452027549574158825829094058509038207331759084319977843276042858638373474688249348076992873314052238942308690507303954406912580102486061607236995039108474881561762586455122590155429569794161571069411135036698761943155931868172074076913102797191291220493392121157");
BigDecimal R = _5.multiply(new BigDecimal("16"),mc);
R = R.subtract(_239.multiply(new BigDecimal("4"),mc),mc);
println(R);

ここで、arctan(1/x)の近似式の収束について考える。

\tan^{-1}\Big(\frac{1}{5}\Big) =
\frac{1}{5}-
\frac{1}{3\times5^{3}}+
\frac{1}{5\times5^{5}}-
\frac{1}{7\times5^{7}}+
\frac{1}{9\times5^{9}}-
\cdots

例えばarctan(1/5)について、20項まで円周率の正しい桁を抜き出す。

項	値
1	0.2
2	0.1973
...
19	0.197395559849880758370049765
20	0.197395559849880758370049765194

20項で約30桁正しく、線形に精度は向上する。
300桁ほしいなら、200項ぐらい計算が必要なはず。

同様にarctan(1/239)について、20項まで円周率の正しい桁を抜き出す。

項	値
1	0.004184
2	0.00418407600
...
19	0.0041840760020747238645382149592854527410480653076319508270196128871817783414228932737826058136
20	0.00418407600207472386453821495928545274104806530763195082701961288718177834142289327378260581362290

20項で約96+2桁正しく、線形に精度は向上する。
300桁ほしいなら、60項ぐらい計算が必要なはず。

同様にarctan(1/49)の場合は、
20項で約70桁正しく、線形に精度は向上する。
300桁ほしいなら、90項ぐらい計算が必要なはず。

arctan(1/X)を計算するコード

無駄が多い気がするが、一応掲載。

import java.math.*;

MathContext mc = new MathContext(100);//必要な桁数以上にする

BigDecimal one = new BigDecimal("1");
BigDecimal mOne = new BigDecimal("-1");
BigDecimal two = new BigDecimal("2");

BigDecimal X = new BigDecimal("5");//必要とする分母の整数

void setup() {
  BigDecimal sum = new BigDecimal("0.0");
  for (int k=0; k<20; k++) {
    sum = sum.add(KOU(k));
    println(sum);
  }
}

BigDecimal KOU(int k) {
  BigDecimal K = new BigDecimal(k);
  BigDecimal bottom;//分母
  BigDecimal top;//分子

  bottom = X.pow(2*k+1);
  bottom = bottom.multiply(K.multiply(two).add(one));
  top = mOne.pow(k);
  return top.divide(bottom, mc);
}

void draw() {
}

arctan型の式のまとめ

できるだけarctan(1/x)のxを大きくするほうが、収束が速い
近似式全体の計算量は、分母が小さい部分に引っ張られると思うので、以下のような式が収束が速いと思う。
それでも、x=49,57,107なので、効果の程はそんなでもない気がするが、未計算。

それにしても、式の美しさと、収束の速さを両立する究極の式ってのは、ないのか？

Stormer

\frac{\pi}{4}=
44\tan^{-1}\Big(\frac{1}{57}\Big)+
7\tan^{-1}\Big(\frac{1}{239}\Big)-
12\tan^{-1}\Big(\frac{1}{682}\Big)+
24\tan^{-1}\Big(\frac{1}{12943}\Big)

高野喜久雄

\frac{\pi}{4}=
12\tan^{-1}\Big(\frac{1}{49}\Big)+
32\tan^{-1}\Big(\frac{1}{57}\Big)-
5\tan^{-1}\Big(\frac{1}{239}\Big)+
12\tan^{-1}\Big(\frac{1}{110443}\Big)

清宮宏

\frac{\pi}{4}=
83\tan^{-1}\Big(\frac{1}{107}\Big)+
17\tan^{-1}\Big(\frac{1}{1710}\Big)-
44\tan^{-1}\Big(\frac{1}{225443}\Big)-
68\tan^{-1}\Big(\frac{1}{2513489}\Big)+
22\tan^{-1}\Big(\frac{1}{42483057}\Big)+
34\tan^{-1}\Big(\frac{1}{7939642926390344818}\Big)

M.Wetherfield

\frac{\pi}{4}=
83\tan^{-1}\Big(\frac{1}{107}\Big)+
17\tan^{-1}\Big(\frac{1}{1710}\Big)-
22\tan^{-1}\Big(\frac{1}{103697}\Big)-
24\tan^{-1}\Big(\frac{1}{2513489}\Big)-
44\tan^{-1}\Big(\frac{1}{18280007883}\Big)+
12\tan^{-1}\Big(\frac{1}{7939642926390344818}\Big)+
22\tan^{-1}\Big(\frac{1}{3054211727257704725384731479018}\Big)

ディリクレ積分

\int_{0}^{\infty}\frac{\sin x}{x}dx=\frac{\pi}{2}

積分を台形式でやってはみたものの、ある程度頑張っても、収束が遅く、3.14以降の精度が芳しくない。

ラマヌジャン1

\frac{1}{\pi}=\frac{\sqrt{8}}{99^2}\sum_{n=0}^\infty \frac{(4n)!}{(4^n n!)^4}\frac{1103+26390n}{99^{4n}}

1項で誤差0.0000001
2項で誤差0.0000000000000001
3項で誤差0.000000000000000000000001
4項で誤差0.00000000000000000000000000000001
８桁ずつ精度が上がる

ラマヌジャン2

\frac{4}{\pi}=\frac{1}{882}\sum_{n=0}^\infty \frac{(-1)^n(4n)!}{(4^n n!)^4}\frac{1123+21460n}{882^{2n}}

1項で誤差0.00001
2項で誤差0.00000000001
６桁ずつ精度が上がる

チェドノフスキー

\frac{1}{\pi}=12\sum_{n=0}^\infty \frac{(-1)^n(6n)!}{(3n)!(n!)^3}\frac{13591409+545140134n}{(640320^3)^{n+\frac{1}{2}}}

1項で誤差0.0000000000001
2項で誤差0.00000000000000000000000001
3項で誤差0.0000000000000000000000000000000000000001
14桁ずつ精度が上がる
この式に対して、計算量の最適化などをしたものが、最速らしい。次の式のほうが収束速そうなんだけど、なぜ？

ガウス-ルジャンドル

下記のpythonコードをBigDecimalを使いprocessingに書き直した。
https://ja.wikipedia.org/wiki/ガウス＝ルジャンドルのアルゴリズム
10万桁を約20秒で計算できる。
100万桁を約2000秒(33分)で計算できる。
とりあえず100桁までは、正しいことを確認した。

1ループで
2.9
２ループ以降、正確な範囲
3.14
3.1415926
3.141592653589793238
3.1415926535897932384626433832795028841971
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862

確かに、ほぼ、倍々で増えていく。すごい。

オリジナルpython版は
10万桁を約13秒で計算できる。
100万桁を約203秒(3.5分)で計算できる。

※javaのBigDecimalはpythonより遅い

super pi 104万桁が10秒とかなので、super piの最適化がすごいとか、桁と計算時間が線形にならないこととか、いろいろ感じるものがある。

calc.py

  import time

  ...

  start = time.time()
  print(picalc())
  end = time.time()
  print(end-start)

calc.pde

import java.math.*;

MathContext mc = new MathContext(100);

BigDecimal a = new BigDecimal(1);
BigDecimal _b = new BigDecimal(2).sqrt(mc);
BigDecimal b = new BigDecimal(1).divide(_b, mc);
BigDecimal _t = new BigDecimal(4);
BigDecimal t = new BigDecimal(1).divide(_t, mc);
BigDecimal p = new BigDecimal(1);
BigDecimal r = new BigDecimal(0);
BigDecimal rn = new BigDecimal(3);
BigDecimal an = new BigDecimal(0);
BigDecimal bn = new BigDecimal(0);
BigDecimal tn = new BigDecimal(0);
BigDecimal pn = new BigDecimal(0);

BigDecimal a_an,a_b,_4t;

while(!(r.equals(rn))){
  r=rn;
  an = a.add(b);
  an = an.divide(new BigDecimal(2), mc);
  bn = a.multiply(b);
  bn = bn.sqrt(mc);
  a_an = a.subtract(an);
  a_an = a_an.multiply(a_an);
  a_an = a_an.multiply(p);
  tn = t.subtract(a_an);
  pn = p.multiply(new BigDecimal(2));
  a_b = a.add(b);
  a_b = a_b.multiply(a_b);
  _4t = t.multiply(new BigDecimal(4));
  rn = a_b.divide(_4t, mc);
  a = an;
  b = bn;
  t = tn;
  p = pn;
}

println(rn);

参考

円周率.jp

近藤 茂
パソコンによる円周率小数点以下５兆桁の計算

円周率の公式集：松元隆二
arctan系の式650万見つけたらしい。

超高速！多倍長整数の計算手法【後編：N! の計算から円周率 100 万桁の挑戦まで】
arctanの近似式解説もある

円周率の公式が多く掲載

ラマヌジャン型公式

ラマヌジャンとチュドノフスキーの比較

計算効率の良いarctan関係式の探索の試み

円周率とarctan型公式

逆三角関数のArcTanの値が欲しい

ライプニッツの公式

定積分

台形式

円周率100桁

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

円周率1000桁

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199

arctan(1/5)100桁

0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441

arctan(1/5)1000桁

0.1973955598498807583700497651947902934475851037878521015176889402410339699782437857326978280372880441126281180736913601044564798867942393557475654952163032700522107470015645015560061286185526633257318692806643896806189528405825931124251613297313993397113233537821796084176648310525473039665725650488878155309384290579311695934192851806364919697519401708560949527368673738508400812367856158009329822514023246675549211026704574378815474839079978985020075223696837961392278354193255722328413846477441352909705465122438302697560518377574220877835853152464749330914587633823112490332030126805100670223312575050942448460267162254894079226140467995062365969282873058287872053603034570766066681374312566267431408992605741703545394046813622411076808100362137599259589071666485145255706799548810043132954668389229036883093278538312521768924841625558744394478718875141137978430670810332362273330875483507230420969155374569324447556530899212193854144645083686720518519002662277446939325408961915860225494493742927

arctan(1/２３９)100桁

0.004184076002074723864538214959285452741048065307631950827019612887181778341422893273782605813622909455

arctan(1/239)1000桁

0.004184076002074723864538214959285452741048065307631950827019612887181778341422893273782605813622909454975450666444863756052458394789311865058922128833092800846271962330773375947634603318473414570331986015454814805992449830211460391253949527607796881558881273397853346518045742548135867464475197910232830977002064652827634653296910481838654356078919591451232220944636862766155208316796426465746551103251034352628244512693556704996844452479043317728393070863140193869519503705864107708558554045223553881423767708365156918252702002293089544950043585440934496440142418724950922838623954553335651171973747020234947597790974695011188854766739795731537093032782113089842583083677190910083909851655104192241678092053268649162667402716844424477315796452027549574158825829094058509038207331759084319977843276042858638373474688249348076992873314052238942308690507303954406912580102486061607236995039108474881561762586455122590155429569794161571069411135036698761943155931868172074076913102797191291220493392121157