円周率の計算
（2015.4.13作成）
迷路のプログラムをいじっていたら、やはり昔からやってきた円周率の計算プログラムが気になって、暫くのめり込んでいます。
現代の計算記録は十兆を超えていて、素人がどうにか出来るレベルではありません。数学も素人なので、数式などを書いて理論から説明できる訳でもありません。まあ、ボロが出ない程度に、今までの経過をかいつまんで書き留めてみたいと思います。
思い返せば会社に入りたての頃（云十年前）、マイコンに興味を持ってマイコンで計算プログラムを動かした事から始まっています。トラ技（1979年頃）にTK80（有名なマイコンボード）で円周率を計算したというレポート記事も掲載され、やはり皆同じ事を考えてるなと思って、自分のプログラムと作り方や性能を比べて楽しんでいました。bitという雑誌にも同様の話題が掲載されていたようです。
マイコン計算で使用されたのは、逆正接（アークタンジェント）の和で表される公式、中でも有名なのはマチン、シムソン（クリンゲンシェルナの公式として知られていましたが、最近に先発見者が見つかったようです）、シュテルマー、ガウス、高野といった方々が発見した公式だったと思います。当時の8ビットマイコン（8080／6800クラス）はそれら公式を使って1000桁を概ね数分くらいで計算するという性能でした。当時のマイコンは除算命令が無かった(乗算命令も無し)ので、除算計算を節約できる公式が有利でした。また多桁乗算（FFTなど）は必要なく、レジスタに収まる程度の値の繰り返し演算で結果が求まるので取っつき易かったといえるでしょう。しかし夫々の公式の計算量は本質的に差がなく、計算桁数を10倍にするには100倍の計算時間がかかるという宿命がありました（（注）現在では違います⇒後述）。後になって登場してきたラマヌジャン、チュドノフスキーの公式も、1項の計算毎に得られる精度はマチンなど逆正接を使った公式よりはるかに大きいですが、計算量の増加率でみると逆正接の公式群と同じ括りにあります。それにこれらは、逆正接の公式と違い、要求桁数精度の平方根や逆数計算が必要なので、一寸敷居が高かったです。
私の場合、マイコン時代（8085＋ROM2KByte＋RAM6KByte＋VRAM64KByte＋CPM80(VRAM52Kを利用した仮想メモリでの52K版)の自作マシン）に十進数計算（数値をBCDコード（'00'～'99'）で処理）で
（ａ）マチンの公式：
（ｂ）シムソンの公式：
π/4 = 8 tan-1(1/10) - tan-1(1/239) - 4 tan-1(1/515)
二進数計算（数値はバイト内で二進数処理0～255）で
（ｃ）シュテルマーの公式：
π/4 = 6 tan-1(1/8) + 2 tan-1(1/57) + tan-1(1/239)
（ｄ）高野の公式：
π/4 = 6 tan-1(1/8) + tan-1(1/32) + 2 tan-1(1/239) - tan-1(1/2943)
夫々をアセンブラでプログラムを作り、少しでも早くなるようにチマチマ修正しては性能を確認していました。最近では、8085/Z80エミュレータ（A～FのFDD相当file 1.44MB×2、720KB×2、640KB×2）を作ってCPM80(64KByte版)をそのままPCで動かして、その上で実行させたりしていました。しかし性能向上も限界は見えており、或るところで放り出しました。
それから少し経てFFTプログラムを作ろうという気力が湧いてきて、解説記事を必死におさらいし、Excel等を使って試し演算するなどして、何とか完成させることができました。
それで漸くラマヌジャン、チュドノフスキーの公式を使った計算をする事が出来ました。この頃には級数公式ではなく、算術幾何平均を使った公式群で、漸化式を1つ計算する毎に精度が2／3／4倍になるような優れものが発表されていて、計算記録の主流となっていました。私もそれらのいくつかを試してみました。
（ａ）2次収束の公式A（サラミンとブレントの式）：
（ｂ）2次収束の公式B（サラミンとブレントの改良版）：
（ｃ）3次収束の公式（ボールウェイン兄弟の式）：
（ｄ）4次収束の公式A（ボールウェイン兄弟の式）：
（ｅ）4次収束の公式B（サラミンとブレントの式の変形）：
結果、計算式が最も簡素な（ｂ）2次収束公式B（サラミンとブレントの改良公式）が、一番性能がでました。4乗根計算が2乗根計算の2倍未満の時間で実行できれば4次収束の公式のほうが早いはずですが、それ以外の付属計算が多く、サラミンとブレント改良公式には及んでいません。
世間の計算記録を見ると、チュドノフスキーの公式を使用してるものがあり、級数計算でなぜ世界記録がだせるのか不思議に思っていましたが、どなたかの円周率計算解説記事を拝見し、級数計算の裏技とでも言うべき考え方を知ることができました。
【級数計算の効率向上（高速化）】
正確には何種類か考え方（作り方）の違いがあるのかも知れませんが、私が理解した方法は至極単純です。
・級数を括弧を使って整理する
・括弧を1つ外して連続する2つの項（偶数番目と奇数番目の項）を1つの項にまとめる。
・全体の項が半分になったら、また最初から上記を繰り返す。
・項が1つになったら、初めて除算を行い、結果（実数値）を求める。
これで計算量の増加率が圧倒的に緩くなります（計算桁数が大きい程、計算回数が少なくなる方式となっています）。
以下の3種類で試して見ました。
（ａ）マチンの公式
（ｂ）ラマヌジャンの公式
迷路のプログラムをいじっていたら、やはり昔からやってきた円周率の計算プログラムが気になって、暫くのめり込んでいます。
現代の計算記録は十兆を超えていて、素人がどうにか出来るレベルではありません。数学も素人なので、数式などを書いて理論から説明できる訳でもありません。まあ、ボロが出ない程度に、今までの経過をかいつまんで書き留めてみたいと思います。
思い返せば会社に入りたての頃（云十年前）、マイコンに興味を持ってマイコンで計算プログラムを動かした事から始まっています。トラ技（1979年頃）にTK80（有名なマイコンボード）で円周率を計算したというレポート記事も掲載され、やはり皆同じ事を考えてるなと思って、自分のプログラムと作り方や性能を比べて楽しんでいました。bitという雑誌にも同様の話題が掲載されていたようです。
マイコン計算で使用されたのは、逆正接（アークタンジェント）の和で表される公式、中でも有名なのはマチン、シムソン（クリンゲンシェルナの公式として知られていましたが、最近に先発見者が見つかったようです）、シュテルマー、ガウス、高野といった方々が発見した公式だったと思います。当時の8ビットマイコン（8080／6800クラス）はそれら公式を使って1000桁を概ね数分くらいで計算するという性能でした。当時のマイコンは除算命令が無かった(乗算命令も無し)ので、除算計算を節約できる公式が有利でした。また多桁乗算（FFTなど）は必要なく、レジスタに収まる程度の値の繰り返し演算で結果が求まるので取っつき易かったといえるでしょう。しかし夫々の公式の計算量は本質的に差がなく、計算桁数を10倍にするには100倍の計算時間がかかるという宿命がありました（（注）現在では違います⇒後述）。後になって登場してきたラマヌジャン、チュドノフスキーの公式も、1項の計算毎に得られる精度はマチンなど逆正接を使った公式よりはるかに大きいですが、計算量の増加率でみると逆正接の公式群と同じ括りにあります。それにこれらは、逆正接の公式と違い、要求桁数精度の平方根や逆数計算が必要なので、一寸敷居が高かったです。
私の場合、マイコン時代（8085＋ROM2KByte＋RAM6KByte＋VRAM64KByte＋CPM80(VRAM52Kを利用した仮想メモリでの52K版)の自作マシン）に十進数計算（数値をBCDコード（'00'～'99'）で処理）で
（ａ）マチンの公式：
π/4 = 4 tan-1(1/5) - tan-1(1/239)
（ｂ）シムソンの公式：
π/4 = 8 tan-1(1/10) - tan-1(1/239) - 4 tan-1(1/515)
二進数計算（数値はバイト内で二進数処理0～255）で
（ｃ）シュテルマーの公式：
π/4 = 6 tan-1(1/8) + 2 tan-1(1/57) + tan-1(1/239)
（ｄ）高野の公式：
π/4 = 6 tan-1(1/8) + tan-1(1/32) + 2 tan-1(1/239) - tan-1(1/2943)
夫々をアセンブラでプログラムを作り、少しでも早くなるようにチマチマ修正しては性能を確認していました。最近では、8085/Z80エミュレータ（A～FのFDD相当file 1.44MB×2、720KB×2、640KB×2）を作ってCPM80(64KByte版)をそのままPCで動かして、その上で実行させたりしていました。しかし性能向上も限界は見えており、或るところで放り出しました。
それから少し経てFFTプログラムを作ろうという気力が湧いてきて、解説記事を必死におさらいし、Excel等を使って試し演算するなどして、何とか完成させることができました。
それで漸くラマヌジャン、チュドノフスキーの公式を使った計算をする事が出来ました。この頃には級数公式ではなく、算術幾何平均を使った公式群で、漸化式を1つ計算する毎に精度が2／3／4倍になるような優れものが発表されていて、計算記録の主流となっていました。私もそれらのいくつかを試してみました。
（ａ）2次収束の公式A（サラミンとブレントの式）：
（ｂ）2次収束の公式B（サラミンとブレントの改良版）：
（ｃ）3次収束の公式（ボールウェイン兄弟の式）：
（ｄ）4次収束の公式A（ボールウェイン兄弟の式）：
（ｅ）4次収束の公式B（サラミンとブレントの式の変形）：
結果、計算式が最も簡素な（ｂ）2次収束公式B（サラミンとブレントの改良公式）が、一番性能がでました。4乗根計算が2乗根計算の2倍未満の時間で実行できれば4次収束の公式のほうが早いはずですが、それ以外の付属計算が多く、サラミンとブレント改良公式には及んでいません。
世間の計算記録を見ると、チュドノフスキーの公式を使用してるものがあり、級数計算でなぜ世界記録がだせるのか不思議に思っていましたが、どなたかの円周率計算解説記事を拝見し、級数計算の裏技とでも言うべき考え方を知ることができました。
【級数計算の効率向上（高速化）】
正確には何種類か考え方（作り方）の違いがあるのかも知れませんが、私が理解した方法は至極単純です。
・級数を括弧を使って整理する
・括弧を1つ外して連続する2つの項（偶数番目と奇数番目の項）を1つの項にまとめる。
・全体の項が半分になったら、また最初から上記を繰り返す。
・項が1つになったら、初めて除算を行い、結果（実数値）を求める。
これで計算量の増加率が圧倒的に緩くなります（計算桁数が大きい程、計算回数が少なくなる方式となっています）。
以下の3種類で試して見ました。
（ａ）マチンの公式
（ｂ）ラマヌジャンの公式
（ｃ）チュドノフスキーの公式
ラマヌジャンの公式を例としてみます。
右辺級数部分は
と整理できます。括弧を1つ除いて級数の隣り合う偶数項と奇数項をまとめて書き直すと、
これを見ると最初に整理した形と同じなので、再度2項ずつをまとめていき、1項だけになるまで繰り返します。最後は
の形になります。ここで初めて割り算（Aの逆数計算）を行います。
実際には、1段ずつを最後までまとめるのではなく、2項の係数を計算できたらすぐにまとめ計算をしてしまい、各項の係数全てをいちいち記憶しないで済むように処理します。
そのロジックを確認するため十進BASICで作成した計算プログラムと、計算結果（2000桁）を例として示します。
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16,2)
DIM b(16,2)
DIM c(16,2)
LET nn=512
LET a(0,0)=1123
LET b(0,0)=(-1)*3*2*1
LET c(0,0)=1
FOR i=1 TO nn-1
LET u=MOD(i,2)
LET a(0,u)=1123+21460*i
LET b(0,u)=(-1)*(4*i+3)*(4*i+2)*(4*i+1)
LET c(0,u)=7056^2*i^3
IF u=1 THEN
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
LET a(j,v)=a(j-1,0)*c(j-1,1)+a(j-1,1)*b(j-1,0)
LET b(j,v)=b(j-1,0)*b(j-1,1)
LET c(j,v)=c(j-1,0)*c(j-1,1)
IF v=0 THEN EXIT FOR
NEXT j
END IF
NEXT i
LET PP=3528*c(j,0)/a(j,0)
PRINT " 3."
LET PP=PP-3
FOR i=0 TO 39
PRINT USING "-----%: ":i*50;
FOR j=0 TO 4
LET pp=pp*10000000000
LET pp1=INT(pp)
LET pp=pp-pp1
PRINT USING "%%%%%%%%%% ":pp1;
NEXT j
PRINT
NEXT i
END
他の級数公式でも同様の方法が適用できて、計算量を削減できます。
c言語に焼直したものを動かして、昔どおり級数を1項ずつ計算する方式と比べてその速さに驚かされました。
しかし、十分検証せずに済ませておいたのですが、最近計算結果を下の桁まで調べたら、精度が全く出ていないことが判明しました。十進BASICでは計算途中の桁落ちがなく計算実行できているので気付きませんでした。c言語で作った方は、途中計算結果を適当な長さで保持して処理を進めるようにしてありますが、その適当な桁数での計算途上で桁落ちなどが起こっているのでしょう、甘かったです。
これで暇つぶしネタが1つ増えてしまいました。
（続く）
（注）上記掲載の公式、式の変形、プログラム等に関して、誤りが含まれている可能性が有りますので、そのままご利用なさらないで下さい。
≪参考文献・記事≫
「円周率の公式と計算法」 大浦拓哉
「円周率の公式集」 編集：松元隆二
「πのａｒｃｔａｎｇｅｎｔ公式集」 柴田昭彦
ラマヌジャンの公式を例としてみます。
右辺級数部分は
と整理できます。括弧を1つ除いて級数の隣り合う偶数項と奇数項をまとめて書き直すと、
これを見ると最初に整理した形と同じなので、再度2項ずつをまとめていき、1項だけになるまで繰り返します。最後は
の形になります。ここで初めて割り算（Aの逆数計算）を行います。
実際には、1段ずつを最後までまとめるのではなく、2項の係数を計算できたらすぐにまとめ計算をしてしまい、各項の係数全てをいちいち記憶しないで済むように処理します。
そのロジックを確認するため十進BASICで作成した計算プログラムと、計算結果（2000桁）を例として示します。
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16,2)
DIM b(16,2)
DIM c(16,2)
LET nn=512
LET a(0,0)=1123
LET b(0,0)=(-1)*3*2*1
LET c(0,0)=1
FOR i=1 TO nn-1
LET u=MOD(i,2)
LET a(0,u)=1123+21460*i
LET b(0,u)=(-1)*(4*i+3)*(4*i+2)*(4*i+1)
LET c(0,u)=7056^2*i^3
IF u=1 THEN
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
LET a(j,v)=a(j-1,0)*c(j-1,1)+a(j-1,1)*b(j-1,0)
LET b(j,v)=b(j-1,0)*b(j-1,1)
LET c(j,v)=c(j-1,0)*c(j-1,1)
IF v=0 THEN EXIT FOR
NEXT j
END IF
NEXT i
LET PP=3528*c(j,0)/a(j,0)
PRINT " 3."
LET PP=PP-3
FOR i=0 TO 39
PRINT USING "-----%: ":i*50;
FOR j=0 TO 4
LET pp=pp*10000000000
LET pp1=INT(pp)
LET pp=pp-pp1
PRINT USING "%%%%%%%%%% ":pp1;
NEXT j
NEXT i
END
3.
0: 1415926535 8979323846 2643383279 5028841971 6939937510
50: 5820974944 5923078164 0628620899 8628034825 3421170679
100: 8214808651 3282306647 0938446095 5058223172 5359408128
150: 4811174502 8410270193 8521105559 6446229489 5493038196
200: 4428810975 6659334461 2847564823 3786783165 2712019091
250: 4564856692 3460348610 4543266482 1339360726 0249141273
300: 7245870066 0631558817 4881520920 9628292540 9171536436
350: 7892590360 0113305305 4882046652 1384146951 9415116094
400: 3305727036 5759591953 0921861173 8193261179 3105118548
450: 0744623799 6274956735 1885752724 8912279381 8301194912
500: 9833673362 4406566430 8602139494 6395224737 1907021798
550: 6094370277 0539217176 2931767523 8467481846 7669405132
600: 0005681271 4526356082 7785771342 7577896091 7363717872
650: 1468440901 2249534301 4654958537 1050792279 6892589235
700: 4201995611 2129021960 8640344181 5981362977 4771309960
750: 5187072113 4999999837 2978049951 0597317328 1609631859
800: 5024459455 3469083026 4252230825 3344685035 2619311881
850: 7101000313 7838752886 5875332083 8142061717 7669147303
900: 5982534904 2875546873 1159562863 8823537875 9375195778
950: 1857780532 1712268066 1300192787 6611195909 2164201989
1000: 3809525720 1065485863 2788659361 5338182796 8230301952
1050: 0353018529 6899577362 2599413891 2497217752 8347913151
1100: 5574857242 4541506959 5082953311 6861727855 8890750983
1150: 8175463746 4939319255 0604009277 0167113900 9848824012
1200: 8583616035 6370766010 4710181942 9555961989 4676783744
1250: 9448255379 7747268471 0404753464 6208046684 2590694912
1300: 9331367702 8989152104 7521620569 6602405803 8150193511
1350: 2533824300 3558764024 7496473263 9141992726 0426992279
1400: 6782354781 6360093417 2164121992 4586315030 2861829745
1450: 5570674983 8505494588 5869269956 9092721079 7509302955
1500: 3211653449 8720275596 0236480665 4991198818 3479775356
1550: 6369807426 5425278625 5181841757 4672890977 7727938000
1600: 8164706001 6145249192 1732172147 7235014144 1973568548
1650: 1613611573 5255213347 5741849468 4385233239 0739414333
1700: 4547762416 8625189835 6948556209 9219222184 2725502542
1750: 5688767179 0494601653 4668049886 2723279178 6085784383
1800: 8279679766 8145410095 3883786360 9506800642 2512520511
1850: 7392984896 0841284886 2694560424 1965285022 2106611863
1900: 0674427862 2039194945 0471237137 8696095636 4371917287
1950: 4677646575 7396241389 0865832645 9958133904 7802759009
c言語に焼直したものを動かして、昔どおり級数を1項ずつ計算する方式と比べてその速さに驚かされました。
しかし、十分検証せずに済ませておいたのですが、最近計算結果を下の桁まで調べたら、精度が全く出ていないことが判明しました。十進BASICでは計算途中の桁落ちがなく計算実行できているので気付きませんでした。c言語で作った方は、途中計算結果を適当な長さで保持して処理を進めるようにしてありますが、その適当な桁数での計算途上で桁落ちなどが起こっているのでしょう、甘かったです。
これで暇つぶしネタが1つ増えてしまいました。
（続く）
≪参考文献・記事≫
「円周率の公式と計算法」 大浦拓哉
「円周率の公式集」 編集：松元隆二
「πのａｒｃｔａｎｇｅｎｔ公式集」 柴田昭彦