円周率の計算（級数形式の計算節約と計算時間の比較）
（2015.11.6作成）
細かい（節約の）話です。
ラマヌジャン、マチンの公式の様な級数形式の計算では、従来の計算方法によれば必要な桁数精度が得られる項で計算を打ち切ればいいのですが、2項ずつまとめる計算方法では、2の累乗個数の項を計算しないと、ピラミッドの積み上げで最後の1項が求まりません。そのため必要精度以上の項を無駄に計算しなくてはなりません。
再度ラマヌジャンの公式を例としてみます。級数は
と整理でき、括弧を1つ除いて級数の隣り合う偶数項と奇数項をまとめると、
2項ずつのまとめ計算を、1項だけになるまで繰り返します。最後に
の形になって結果が求まります。
ここで、必要精度が求まる項より後の項については、An=0、Bn=1、Cn=1として計算します。まとめ計算の相手が0や1の場合、実際の加算や乗算を省略でき、計算時間を短縮できます。
これを十進BASICで作ると次のようになります（偶数項、奇数項ともに配列に記憶する方法で作成してあります）。
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16,2)
DIM b(16,2)
DIM c(16,2)
LET nd=5000
LET n0=INT(nd/5.8)
LET nn=2^INT(LOG2(n0)+1)
PRINT "LOOP ";n0;"/";nn
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)
IF i<=n0 THEN
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
ELSE
LET a(0,u)=0
LET b(0,u)=1
LET c(0,u)=1
END IF
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
LET iz=INT(nd/50)-1
FOR i=0 TO iz
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
5000桁を求めるためには概ね862項(1000桁を求めるためには172項)までの計算で済むので、それ以降は節約計算を行います。
ただし十進BASICの言語処理では、これが計算時間の節約になっているかどうかは分かりません。
これをチュドノフスキーの公式を使用した計算（2^281*125 + 1を法とする整数環ＦＦＴとKARATSUBAの併用）に組み込んで、約8000万桁の計算に約44分かかりましたが、以前の計算では約1時間かかりましたので、25%計算が節約できたことになります。
使用PC： PCショップオリジナル（WinXP(32bit)、CPUはAMD Athlon 64X2 Dual Core Proc 5400+ CLK 2.8GHz）
ついでに、他のプログラムと比較してみました（自作版は全て、2^281*125 + 1を法とする整数環ＦＦＴとKARATSUBAの併用による計算）。
・自作 チュドノフスキーの公式の計算 2000万桁： 約500秒
・自作 サラミンとブレントの公式（2次収束 改良版）の計算 2000万桁： 約1510秒
・自作 サラミンとブレントの公式（4次収束 変形版）の計算 2000万桁： 約1580秒
・Super π for Windows(Ver 1.1) （©Kanada Labo. University of Tokyo）1677万桁： 約850秒
・y-cruncher v0.5.5 Build 9180 (fix 2) ( www.numberworld.org )
© 2008-2011 Alexander J. Yee ( a-yee@u.northwestern.edu ) 2500万桁： 約120秒
素人の作りとしてまずまずではないでしょうか。
まだやれる事は残っていますが、ここで一旦区切りをつけておきます。強烈な手法が学べたら、改めて挑戦することにします。
（注）上記掲載の公式、式の変形、プログラム等に関して、誤りが含まれている可能性が有りますので、そのままご利用なさらないで下さい。
≪参考文献・記事≫
「円周率の公式と計算法」 大浦拓哉
「円周率の公式集」 編集：松元隆二
細かい（節約の）話です。
ラマヌジャン、マチンの公式の様な級数形式の計算では、従来の計算方法によれば必要な桁数精度が得られる項で計算を打ち切ればいいのですが、2項ずつまとめる計算方法では、2の累乗個数の項を計算しないと、ピラミッドの積み上げで最後の1項が求まりません。そのため必要精度以上の項を無駄に計算しなくてはなりません。
再度ラマヌジャンの公式を例としてみます。級数は
と整理でき、括弧を1つ除いて級数の隣り合う偶数項と奇数項をまとめると、
2項ずつのまとめ計算を、1項だけになるまで繰り返します。最後に
の形になって結果が求まります。
ここで、必要精度が求まる項より後の項については、An=0、Bn=1、Cn=1として計算します。まとめ計算の相手が0や1の場合、実際の加算や乗算を省略でき、計算時間を短縮できます。
これを十進BASICで作ると次のようになります（偶数項、奇数項ともに配列に記憶する方法で作成してあります）。
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16,2)
DIM b(16,2)
DIM c(16,2)
LET nd=5000
LET n0=INT(nd/5.8)
LET nn=2^INT(LOG2(n0)+1)
PRINT "LOOP ";n0;"/";nn
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)
IF i<=n0 THEN
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
ELSE
LET a(0,u)=0
LET b(0,u)=1
LET c(0,u)=1
END IF
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
LET iz=INT(nd/50)-1
FOR i=0 TO iz
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
LOOP 862 / 1024
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
2000: 9465764078 9512694683 9835259570 9825822620 5224894077
2050: 2671947826 8482601476 9909026401 3639443745 5305068203
2100: 4962524517 4939965143 1429809190 6592509372 2169646151
2150: 5709858387 4105978859 5977297549 8930161753 9284681382
2200: 6868386894 2774155991 8559252459 5395943104 9972524680
2250: 8459872736 4469584865 3836736222 6260991246 0805124388
2300: 4390451244 1365497627 8079771569 1435997700 1296160894
2350: 4169486855 5848406353 4220722258 2848864815 8456028506
2400: 0168427394 5226746767 8895252138 5225499546 6672782398
2450: 6456596116 3548862305 7745649803 5593634568 1743241125
2500: 1507606947 9451096596 0940252288 7971089314 5669136867
2550: 2287489405 6010150330 8617928680 9208747609 1782493858
2600: 9009714909 6759852613 6554978189 3129784821 6829989487
2650: 2265880485 7564014270 4775551323 7964145152 3746234364
2700: 5428584447 9526586782 1051141354 7357395231 1342716610
2750: 2135969536 2314429524 8493718711 0145765403 5902799344
2800: 0374200731 0578539062 1983874478 0847848968 3321445713
2850: 8687519435 0643021845 3191048481 0053706146 8067491927
2900: 8191197939 9520614196 6342875444 0643745123 7181921799
2950: 9839101591 9561814675 1426912397 4894090718 6494231961
3000: 5679452080 9514655022 5231603881 9301420937 6213785595
3050: 6638937787 0830390697 9207734672 2182562599 6615014215
3100: 0306803844 7734549202 6054146659 2520149744 2850732518
3150: 6660021324 3408819071 0486331734 6496514539 0579626856
3200: 1005508106 6587969981 6357473638 4052571459 1028970641
3250: 4011097120 6280439039 7595156771 5770042033 7869936007
3300: 2305587631 7635942187 3125147120 5329281918 2618612586
3350: 7321579198 4148488291 6447060957 5270695722 0917567116
3400: 7229109816 9091528017 3506712748 5832228718 3520935396
3450: 5725121083 5791513698 8209144421 0067510334 6711031412
3500: 6711136990 8658516398 3150197016 5151168517 1437657618
3550: 3515565088 4909989859 9823873455 2833163550 7647918535
3600: 8932261854 8963213293 3089857064 2046752590 7091548141
3650: 6549859461 6371802709 8199430992 4488957571 2828905923
3700: 2332609729 9712084433 5732654893 8239119325 9746366730
3750: 5836041428 1388303203 8249037589 8524374417 0291327656
3800: 1809377344 4030707469 2112019130 2033038019 7621101100
3850: 4492932151 6084244485 9637669838 9522868478 3123552658
3900: 2131449576 8572624334 4189303968 6426243410 7732269780
3950: 2807318915 4411010446 8232527162 0105265227 2111660396
4000: 6655730925 4711055785 3763466820 6531098965 2691862056
4050: 4769312570 5863566201 8558100729 3606598764 8611791045
4100: 3348850346 1136576867 5324944166 8039626579 7877185560
4150: 8455296541 2665408530 6143444318 5867697514 5661406800
4200: 7002378776 5913440171 2749470420 5622305389 9456131407
4250: 1127000407 8547332699 3908145466 4645880797 2708266830
4300: 6343285878 5698305235 8089330657 5740679545 7163775254
4350: 2021149557 6158140025 0126228594 1302164715 5097925923
4400: 0990796547 3761255176 5675135751 7829666454 7791745011
4450: 2996148903 0463994713 2962107340 4375189573 5961458901
4500: 9389713111 7904297828 5647503203 1986915140 2870808599
4550: 0480109412 1472213179 4764777262 2414254854 5403321571
4600: 8530614228 8137585043 0633217518 2979866223 7172159160
4650: 7716692547 4873898665 4949450114 6540628433 6639379003
4700: 9769265672 1463853067 3609657120 9180763832 7166416274
4750: 8888007869 2560290228 4721040317 2118608204 1900042296
4800: 6171196377 9213375751 1495950156 6049631862 9472654736
4850: 4252308177 0367515906 7350235072 8354056704 0386743513
4900: 6222247715 8915049530 9844489333 0963408780 7693259939
4950: 7805419341 4473774418 4263129860 8099888687 4132604721
ただし十進BASICの言語処理では、これが計算時間の節約になっているかどうかは分かりません。
これをチュドノフスキーの公式を使用した計算（2^281*125 + 1を法とする整数環ＦＦＴとKARATSUBAの併用）に組み込んで、約8000万桁の計算に約44分かかりましたが、以前の計算では約1時間かかりましたので、25%計算が節約できたことになります。
使用PC： PCショップオリジナル（WinXP(32bit)、CPUはAMD Athlon 64X2 Dual Core Proc 5400+ CLK 2.8GHz）
ついでに、他のプログラムと比較してみました（自作版は全て、2^281*125 + 1を法とする整数環ＦＦＴとKARATSUBAの併用による計算）。
・自作 チュドノフスキーの公式の計算 2000万桁： 約500秒
・自作 サラミンとブレントの公式（2次収束 改良版）の計算 2000万桁： 約1510秒
・自作 サラミンとブレントの公式（4次収束 変形版）の計算 2000万桁： 約1580秒
・Super π for Windows(Ver 1.1) （©Kanada Labo. University of Tokyo）1677万桁： 約850秒
・y-cruncher v0.5.5 Build 9180 (fix 2) ( www.numberworld.org )
© 2008-2011 Alexander J. Yee ( a-yee@u.northwestern.edu ) 2500万桁： 約120秒
素人の作りとしてまずまずではないでしょうか。
まだやれる事は残っていますが、ここで一旦区切りをつけておきます。強烈な手法が学べたら、改めて挑戦することにします。
（注）上記掲載の公式、式の変形、プログラム等に関して、誤りが含まれている可能性が有りますので、そのままご利用なさらないで下さい。
≪参考文献・記事≫
「円周率の公式と計算法」 大浦拓哉
「円周率の公式集」 編集：松元隆二