円周率の計算（続き： マチンの公式での高速化）
（2015.4.15作成）
マチンの公式（逆正接）は、級数の形がラマヌジャン系の形と異なるので、整理の形も一寸違います。
蛇足になりますが、自分用のメモなので、一応記録しておきます。
マチンの公式は、
右辺の逆正接関数の級数展開に1／5（または1／239）を掛けて整理すると、
級数の隣り合う偶数項と奇数項をまとめて書き直すと、
これで最初の整理した形と同じになるので、加減算の符号だけ考慮して再度２項ずつをまとめていき、1項だけになるまで繰り返します。最後は
の形になります。ここで初めて割り算（B*Cの逆数計算）を行います。
実際の処理はラマヌジャンの時と同様に行います。
下記はこのロジックを確認するための十進BASICで作成した計算プログラムです、計算結果（1000桁）は同じなので省略します。1段目の係数計算(代入式)は、省略して2段目に式として組み込めますが、上記の変形と対応させるために残してあります。
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16,2)
DIM b(16,2)
DIM c(16,2)
LET nn1=1024
LET nn2=256
FOR i=0 TO nn1-1
LET u=MOD(i,2)
LET b(0,u)=2*i+1
LET c(0,u)=5*5
IF u=1 THEN
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
IF j=1 THEN
LET a(j,v)=b(j-1,1)*c(j-1,1)-b(j-1,0)
ELSE
LET a(j,v)=a(j-1,0)*b(j-1,1)*c(j-1,1)+a(j-1,1)*b(j-1,0)
END IF
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=80*a(j,0)/b(j,0)/c(j,0)
FOR i=0 TO nn2-1
LET u=MOD(i,2)
LET b(0,u)=2*i+1
LET c(0,u)=239*239
IF u=1 THEN
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
IF j=1 THEN
LET a(j,v)=b(j-1,1)*c(j-1,1)-b(j-1,0)
ELSE
LET a(j,v)=a(j-1,0)*b(j-1,1)*c(j-1,1)+a(j-1,1)*b(j-1,0)
END IF
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 QQ=956*a(j,0)/b(j,0)/c(j,0)
LET PP=PP-QQ
PRINT " 3."
LET PP=PP-3
FOR i=0 TO 19
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
（注）上記掲載の公式、式の変形、プログラム等に関して、誤りが含まれている可能性が有りますので、そのままご利用なさらないで下さい。
≪参考文献・記事≫
「円周率の公式と計算法」 大浦拓哉
「円周率の公式集」 編集：松元隆二
「πのａｒｃｔａｎｇｅｎｔ公式集」 柴田昭彦
マチンの公式（逆正接）は、級数の形がラマヌジャン系の形と異なるので、整理の形も一寸違います。
蛇足になりますが、自分用のメモなので、一応記録しておきます。
マチンの公式は、
右辺の逆正接関数の級数展開に1／5（または1／239）を掛けて整理すると、
級数の隣り合う偶数項と奇数項をまとめて書き直すと、
これで最初の整理した形と同じになるので、加減算の符号だけ考慮して再度２項ずつをまとめていき、1項だけになるまで繰り返します。最後は
の形になります。ここで初めて割り算（B*Cの逆数計算）を行います。
実際の処理はラマヌジャンの時と同様に行います。
下記はこのロジックを確認するための十進BASICで作成した計算プログラムです、計算結果（1000桁）は同じなので省略します。1段目の係数計算(代入式)は、省略して2段目に式として組み込めますが、上記の変形と対応させるために残してあります。
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16,2)
DIM b(16,2)
DIM c(16,2)
LET nn1=1024
LET nn2=256
FOR i=0 TO nn1-1
LET u=MOD(i,2)
LET b(0,u)=2*i+1
LET c(0,u)=5*5
IF u=1 THEN
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
IF j=1 THEN
LET a(j,v)=b(j-1,1)*c(j-1,1)-b(j-1,0)
ELSE
LET a(j,v)=a(j-1,0)*b(j-1,1)*c(j-1,1)+a(j-1,1)*b(j-1,0)
END IF
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=80*a(j,0)/b(j,0)/c(j,0)
FOR i=0 TO nn2-1
LET u=MOD(i,2)
LET b(0,u)=2*i+1
LET c(0,u)=239*239
IF u=1 THEN
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
IF j=1 THEN
LET a(j,v)=b(j-1,1)*c(j-1,1)-b(j-1,0)
ELSE
LET a(j,v)=a(j-1,0)*b(j-1,1)*c(j-1,1)+a(j-1,1)*b(j-1,0)
END IF
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 QQ=956*a(j,0)/b(j,0)/c(j,0)
LET PP=PP-QQ
PRINT " 3."
LET PP=PP-3
FOR i=0 TO 19
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
（注）上記掲載の公式、式の変形、プログラム等に関して、誤りが含まれている可能性が有りますので、そのままご利用なさらないで下さい。
≪参考文献・記事≫
「円周率の公式と計算法」 大浦拓哉
「円周率の公式集」 編集：松元隆二
「πのａｒｃｔａｎｇｅｎｔ公式集」 柴田昭彦