円周率の計算（続き： ラマヌジャンとＢＢＰの公式、並列処理）
（2015.6.16作成）
級数公式での計算用メモリを少し削減してみました（大した話ではありません、蛇足の範囲です）。
奇数項は計算した後、直ぐに上位でのまとめ計算をして用済みになるので、各段の保持変数は偶数項だけにしてメモリを節約しました。
変更したラマヌジャンの公式の、十進BASICの計算プログラムです。（計算結果は省略）
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16)
DIM b(16)
DIM c(16)
LET nn=256
LET a(0)=1123
LET b(0)=(-1)*3*2*1
LET c(0)=1
FOR i=1 TO nn-1
LET u=MOD(i,2)
IF u=0 THEN
LET a(0)=1123+21460*i
LET b(0)=(-1)*(4*i+3)*(4*i+2)*(4*i+1)
LET c(0)=7056^2*i^3
ELSE
LET a1=1123+21460*i
LET b1=(-1)*(4*i+3)*(4*i+2)*(4*i+1)
LET c1=7056^2*i^3
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
IF v=0 THEN
LET a(j)=a(j-1)*c1+a1*b(j-1)
LET b(j)=b(j-1)*b1
LET c(j)=c(j-1)*c1
ELSE
LET a1=a(j-1)*c1+a1*b(j-1)
LET b1=b(j-1)*b1
LET c1=c(j-1)*c1
END IF
IF v=0 THEN EXIT FOR
NEXT j
END IF
NEXT i
LET PP=3528*c(j)/a(j)
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億桁超の計算メモリはとれませんでした（残念）。
速度向上については、累乗根計算処理を少しいじってみました。4乗根の計算が少し早くなり、サラミンとブレントの2次収束（改良）の公式と4次収束の公式の速度がトントンになりました（でも、性能はチュドノフスキーの公式の半分以下）。
一応1億桁超という事で、速度は劣りますがサラミンとブレントの4次収束の公式で約1.6億桁計算してみました。さすがに時間かかりました（約4時間50分）。確認できないと計算が無駄骨になってしまうので、チュドノフスキーの公式での計算（16進で約6700万桁）結果と16進数で比較したら、最後の100ビット位を除いて一致しました。
16進数での任意桁を計算できるBBPの公式というのがあるので、計算プログラムを作ってみました。BBP（16進32桁を計算）との比較でも、16進67108736桁（8000万桁精度相当）と、134217600桁（1億6000万桁精度相当）で一致しました。
これで一応1億桁超まで計算できたことになりました。
BBP公式の確認の為に十進BASICで作成したプログラム（16進8桁表示計算の例）です（参考記事どおりの作りです）。
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE EXTERNAL FUNCTION cmod
LET nn=1024
LET p1=0
FOR i=0 TO nn
LET n1=nn-i
LET p1=p1+cmod(4,8*i+1,n1)/(8*i+1)
LET p1=p1-cmod(2,8*i+4,n1)/(8*i+4)
LET p1=p1-cmod(1,8*i+5,n1)/(8*i+5)
LET p1=p1-cmod(1,8*i+6,n1)/(8*i+6)
NEXT i
FOR j=1 TO 16
LET p1=p1+(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6))/16^j
LET i=i+1
NEXT j
LET p1=p1-INT(p1)
LET p1=INT(p1*16^8)
PRINT BSTR$(p1,16)
END
EXTERNAL FUNCTION cmod(h,k,n)
OPTION ARITHMETIC DECIMAL_HIGH
LET m=MOD(16,k)
FOR j=0 TO 30
IF MOD(n,2)=1 THEN LET h=MOD(h*m,k)
LET n=INT(n/2)
IF n=0 THEN EXIT FOR
LET m=MOD(m*m,k)
NEXT j
LET cmod=h
END FUNCTION
昔マイコンボード作りに没頭していた頃、丁度ＮＥＣ製スパコンSX-2で1億桁を計算したというニュースを聞いたのを思い出しました。それから程なく日立のスパコンで記録が更新されたとか．．．。あの当時のスパコンと今の5万程度のノートPCが同じくらいの性能とは、感慨深いものがあります。
他の対策としては、64ビットモードで書き直してプログラムを動かせばそれなりに早くなりそうですが、何か高性能ＣＰＵに取替えたのと一緒のような感じ（早くなるのが当たり前的）で、イマイチ乗り気になりません。
級数系公式の計算では各項を単独に計算できるので並列計算向きです。効率化したピラミッドの積み上げ計算にしても同様で、適当に分割することで簡単に並列処理する事ができます。分割した計算をスレッド（または別プログラムや別のPC）に割り当てれば、最後のまとめ計算と逆数計算などを除いて単純に並列分だけ早くなるでしょう。
チュドノフスキーの公式（8000万桁計算）で試してみました。ループカウンタを前半後半に2分割して夫々を別プログラムで動かします。
前半分（1/2）をAP1、後半分（2/2）をAP2、最終計算をAP3とします。中間結果はファイルで引き継ぎます。
AP1を単独で動かしたら1735秒、AP1とAP2を同時実行させたら、AP1が1810秒、AP2が1722秒、そして両方が終了したところでまとめ計算のAP3を実行して、465秒かかりました。全てを直列で動かすと約1時間程かかりますが、2分割並列で約38分に短縮できました。CPUコアが4以上なら25分位に短縮できるでしょう。
算術幾何平均系の公式は漸化式なので、並列計算は素人が安直には思いつけないので、級数公式より不利ですね。
（注）上記掲載の内容に関して、誤りが含まれている可能性が有りますので、そのままご利用なさらないで下さい。
≪参考文献・記事≫
「円周率の公式集」 編集：松元隆二
円周率.jp （ウィキペディア）
級数公式での計算用メモリを少し削減してみました（大した話ではありません、蛇足の範囲です）。
奇数項は計算した後、直ぐに上位でのまとめ計算をして用済みになるので、各段の保持変数は偶数項だけにしてメモリを節約しました。
変更したラマヌジャンの公式の、十進BASICの計算プログラムです。（計算結果は省略）
OPTION BASE 0
OPTION ARITHMETIC RATIONAL
DIM a(16)
DIM b(16)
DIM c(16)
LET nn=256
LET a(0)=1123
LET b(0)=(-1)*3*2*1
LET c(0)=1
FOR i=1 TO nn-1
LET u=MOD(i,2)
IF u=0 THEN
LET a(0)=1123+21460*i
LET b(0)=(-1)*(4*i+3)*(4*i+2)*(4*i+1)
LET c(0)=7056^2*i^3
ELSE
LET a1=1123+21460*i
LET b1=(-1)*(4*i+3)*(4*i+2)*(4*i+1)
LET c1=7056^2*i^3
FOR j=1 TO 15
LET v=MOD(INT(i/2^j),2)
IF v=0 THEN
LET a(j)=a(j-1)*c1+a1*b(j-1)
LET b(j)=b(j-1)*b1
LET c(j)=c(j-1)*c1
ELSE
LET a1=a(j-1)*c1+a1*b(j-1)
LET b1=b(j-1)*b1
LET c1=c(j-1)*c1
END IF
IF v=0 THEN EXIT FOR
NEXT j
END IF
NEXT i
LET PP=3528*c(j)/a(j)
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
この方法をチュドノフスキーの公式計算に適用して動かしましたが、まだ我が家のＰＣでは1億桁超の計算メモリはとれませんでした（残念）。
速度向上については、累乗根計算処理を少しいじってみました。4乗根の計算が少し早くなり、サラミンとブレントの2次収束（改良）の公式と4次収束の公式の速度がトントンになりました（でも、性能はチュドノフスキーの公式の半分以下）。
一応1億桁超という事で、速度は劣りますがサラミンとブレントの4次収束の公式で約1.6億桁計算してみました。さすがに時間かかりました（約4時間50分）。確認できないと計算が無駄骨になってしまうので、チュドノフスキーの公式での計算（16進で約6700万桁）結果と16進数で比較したら、最後の100ビット位を除いて一致しました。
16進数での任意桁を計算できるBBPの公式というのがあるので、計算プログラムを作ってみました。BBP（16進32桁を計算）との比較でも、16進67108736桁（8000万桁精度相当）と、134217600桁（1億6000万桁精度相当）で一致しました。
これで一応1億桁超まで計算できたことになりました。
BBP公式の確認の為に十進BASICで作成したプログラム（16進8桁表示計算の例）です（参考記事どおりの作りです）。
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE EXTERNAL FUNCTION cmod
LET nn=1024
LET p1=0
FOR i=0 TO nn
LET n1=nn-i
LET p1=p1+cmod(4,8*i+1,n1)/(8*i+1)
LET p1=p1-cmod(2,8*i+4,n1)/(8*i+4)
LET p1=p1-cmod(1,8*i+5,n1)/(8*i+5)
LET p1=p1-cmod(1,8*i+6,n1)/(8*i+6)
NEXT i
FOR j=1 TO 16
LET p1=p1+(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6))/16^j
LET i=i+1
NEXT j
LET p1=p1-INT(p1)
LET p1=INT(p1*16^8)
PRINT BSTR$(p1,16)
END
EXTERNAL FUNCTION cmod(h,k,n)
OPTION ARITHMETIC DECIMAL_HIGH
LET m=MOD(16,k)
FOR j=0 TO 30
IF MOD(n,2)=1 THEN LET h=MOD(h*m,k)
LET n=INT(n/2)
IF n=0 THEN EXIT FOR
LET m=MOD(m*m,k)
NEXT j
LET cmod=h
END FUNCTION
昔マイコンボード作りに没頭していた頃、丁度ＮＥＣ製スパコンSX-2で1億桁を計算したというニュースを聞いたのを思い出しました。それから程なく日立のスパコンで記録が更新されたとか．．．。あの当時のスパコンと今の5万程度のノートPCが同じくらいの性能とは、感慨深いものがあります。
他の対策としては、64ビットモードで書き直してプログラムを動かせばそれなりに早くなりそうですが、何か高性能ＣＰＵに取替えたのと一緒のような感じ（早くなるのが当たり前的）で、イマイチ乗り気になりません。
級数系公式の計算では各項を単独に計算できるので並列計算向きです。効率化したピラミッドの積み上げ計算にしても同様で、適当に分割することで簡単に並列処理する事ができます。分割した計算をスレッド（または別プログラムや別のPC）に割り当てれば、最後のまとめ計算と逆数計算などを除いて単純に並列分だけ早くなるでしょう。
チュドノフスキーの公式（8000万桁計算）で試してみました。ループカウンタを前半後半に2分割して夫々を別プログラムで動かします。
前半分（1/2）をAP1、後半分（2/2）をAP2、最終計算をAP3とします。中間結果はファイルで引き継ぎます。
AP1を単独で動かしたら1735秒、AP1とAP2を同時実行させたら、AP1が1810秒、AP2が1722秒、そして両方が終了したところでまとめ計算のAP3を実行して、465秒かかりました。全てを直列で動かすと約1時間程かかりますが、2分割並列で約38分に短縮できました。CPUコアが4以上なら25分位に短縮できるでしょう。
算術幾何平均系の公式は漸化式なので、並列計算は素人が安直には思いつけないので、級数公式より不利ですね。
（注）上記掲載の内容に関して、誤りが含まれている可能性が有りますので、そのままご利用なさらないで下さい。
≪参考文献・記事≫
「円周率の公式集」 編集：松元隆二
円周率.jp （ウィキペディア）