円周率の計算（その４： 乗算方法など）
（2015.5.27作成）
整数環ＦＦＴ、私のPCでは億の桁長を計算できないので、ビット数を削って若干速度を改善しました。
例８．P=2151*315 + 1 (16進表示: 0x9D800000, 0, 0, 0, 1)
(10進表示: 899166046404754725066720160753181935921130373121)
ω=97 (w=2150)
1桁の単位=64ビット
例９．P=2281*125 + 1 (16進表示: 0xFA000000, 0, 0, 0, 0, 0, 0, 0, 1)
(10進表示: 485667223056432267729865476705879726660601709763034880312953102434726071301302124544001)
ω=167 (w=2281)
1桁の単位=128ビット
前回の例６（P=2189*3 + 1、1桁の単位=64ビット)は192ビット乗算を行いますが、例８では160ビットで済むのでその分若干早くなりました。
同様に例７（P=2316*13 + 1、1桁の単位=128ビット)は320ビット乗算ですが、例９では288ビットで済みます。
例８、９で、64ビット浮動小数点計算の複素ＦＦＴと大体同じ位の速度になりました。
現状ではここ辺りが限界でしょうか。
円周率の計算時間は、サラミンとブレントの改良公式はチュドノフスキーの公式での２倍強かかります。
チュドノフスキーの公式の計算では、２５０万桁が約５０秒、．．．、４０００万桁が約２５分、８０００万桁が１時間弱かかりました。これにはＫＡＲＡＴＳＵＢＡ方式も併用（小さい桁の乗算で）しています。
チュドノフスキーの公式では殆どの計算が乗算なので、乗算速度が上がらないと性能が向上しません。しかし算術幾何平均を使った公式群では冪乗根の計算が主体なので、乗算速度が同じでもニュートン法の計算手順のなかでの改善で、見直せる事があるかもしれません（でも２倍の差は埋まらないでしょう、多分）。
しかしその前に計算結果の出力（２進１０進変換）処理の高速化が必要になりました。こちらは全然性能など考えていませんでしたが、桁数が上がるにつれて２進１０進変換時間の方が長くかかる様になってしまいました。なので前記桁数は計算確定の時間ではありませんので、悪しからず。
サラミンとブレントの改良公式で１０００万桁、チュドノフスキーの公式で５００万桁までは、なんとか出力して確認できました。
他に、ワラを掴む的な準備（下調べ）として、Toom-Cook法（４分割）の計算を確認してみました。参考記事では行列表記でさらりと説明されていますが、行列（式）は得意ではないので普通の式に置き換えて、十進ＢＡＳＩＣで計算してみました。すると計算結果が異なる場合があり、正誤のいくつかを見比べて、どうも係数が違っている（らしい？）事が分かりました。正解はｘｘ元連立方程式を解かないと出てこないらしいので、それは難しそうだからスキップして、正しい計算結果となるような係数を逆算で推定してみました（”７２”→”１８００”）。
性能は度外視して、理論通り？と思われる式で計算した結果が次のとおりです。KARATSUBA方式の倍に分割しただけなのに、なんでこんなに式が多くなるんだか．．．。
まだやれそうな（やらなきゃならない）事があるので、取敢えず良し々々。
≪プラグラム例と結果≫
（１６桁×１６桁を４桁ずつに分解して計算しています）
OPTION ARITHMETIC RATIONAL
LET dn=4
LET d0=10^4
LET xx=9988776655443322
LET yy=9999222288881111
PRINT "xx =";xx,"yy =";yy
PRINT "非分割乗算 xx * yy =";
PRINT USING REPEAT$("-",8*dn-1) & "%":xx*yy
PRINT
LET x0=MOD(xx,d0)
LET x1=MOD(INT(xx/d0),d0)
LET x2=MOD(INT(xx/d0^2),d0)
LET x3=MOD(INT(xx/d0^3),d0)
LET y0=MOD(yy,d0)
LET y1=MOD(INT(yy/d0),d0)
LET y2=MOD(INT(yy/d0^2),d0)
LET y3=MOD(INT(yy/d0^3),d0)
PRINT "xx=";x3;x2;x1;x0,"yy=";y3;y2;y1;y0
LET u6= -x3 +2*x2 -4*x1 +8*x0
LET u5=-8*x3 +4*x2 -2*x1 +x0
LET u4= -x3 +x2 -x1 +x0
LET u3= x0
LET u2= x3 +x2 +x1 +x0
LET u1= 8*x3 +4*x2 +2*x1 +x0
LET u0= x3
LET v6= -y3 +2*y2 -4*y1 +8*y0
LET v5=-8*y3 +4*y2 -2*y1 +y0
LET v4= -y3 +y2 -y1 +y0
LET v3= y0
LET v2= y3 +y2 +y1 +y0
LET v1= 8*y3 +4*y2 +2*y1 +y0
LET v0= y3
PRINT "u=";u6;u5;u4;u3;u2;u1;u0
PRINT "v=";v6;v5;v4;v3;v2;v1;v0
LET uv6=u6*v6
LET uv5=u5*v5
LET uv4=u4*v4
LET uv3=u3*v3
LET uv2=u2*v2
LET uv1=u1*v1
LET uv0=u0*v0
PRINT "uv=";uv6;uv5;uv4;uv3;uv2;uv1;uv0
LET w6= 360*uv0
LET w5= -4*uv6-10*uv5+120*uv4+180*uv3 -40*uv2 +6*uv1 +180*uv0
LET w4= 15*uv5 -60*uv4 +90*uv3 -60*uv2 +15*uv1-1800*uv0
LET w3= 20*uv6+20*uv5-540*uv4-900*uv3+140*uv2 -900*uv0
LET w2= -15*uv5+240*uv4-450*uv3+240*uv2 -15*uv1+1440*uv0
LET w1=-16*uv6-10*uv5+240*uv4+720*uv3 +80*uv2 -6*uv1 +720*uv0
LET w0= 360*uv3
LET w6=w6/360
LET w5=w5/360
LET w4=w4/360
LET w3=w3/360
LET w2=w2/360
LET w1=w1/360
LET w0=w0/360
PRINT "w=";w6;w5;w4;w3;w2;w1;w0
LET w1=w1+INT(w0/d0)
LET w0=MOD(w0,d0)
LET w2=w2+INT(w1/d0)
LET w1=MOD(w1,d0)
LET w3=w3+INT(w2/d0)
LET w2=MOD(w2,d0)
LET w4=w4+INT(w3/d0)
LET w3=MOD(w3,d0)
LET w5=w5+INT(w4/d0)
LET w4=MOD(w4,d0)
LET w6=w6+INT(w5/d0)
LET w5=MOD(w5,d0)
LET w7=w7+INT(w6/d0)
LET w6=MOD(w6,d0)
LET fmt1$=REPEAT$("-",dn-1) & "%"
LET fmt2$=REPEAT$("%",dn)
PRINT "Toom-Cook法(４分割) zz = ";
PRINT USING fmt1$:w7;
PRINT USING fmt2$:w6;
PRINT USING fmt2$:w5;
PRINT USING fmt2$:w4;
PRINT USING fmt2$:w3;
PRINT USING fmt2$:w2;
PRINT USING fmt2$:w1;
PRINT USING fmt2$:w0;
END
xx = 9988776655443322 yy = 9999222288881111
非分割乗算 xx * yy =99879998171764182850815056890742
xx= 9988 7766 5544 3322 yy= 9999 2222 8888 1111
u= 9944 -56606 -4444 3322 26620 125378 9988
v=-32219 -87769 -15554 1111 22220 107767 9999
uv=-320385736 4968252014 69121976 3690742 591496400 13511610926 99870012
w= 99870012 99845570 161463852 125656322 65284582 35685320 3690742
Toom-Cook法(４分割) zz = 99879998171764182850815056890742
（注）本文は、数学素人の個人的な記録（単なるメモ）なので、誤りがある事などをご了解ください。
≪参考文献・記事≫
計算代数（加古富志雄H26.4.24）
円周率．ｊｐ （ウィキペディア）
整数環ＦＦＴ、私のPCでは億の桁長を計算できないので、ビット数を削って若干速度を改善しました。
例８．P=2151*315 + 1 (16進表示: 0x9D800000, 0, 0, 0, 1)
(10進表示: 899166046404754725066720160753181935921130373121)
ω=97 (w=2150)
1桁の単位=64ビット
例９．P=2281*125 + 1 (16進表示: 0xFA000000, 0, 0, 0, 0, 0, 0, 0, 1)
(10進表示: 485667223056432267729865476705879726660601709763034880312953102434726071301302124544001)
ω=167 (w=2281)
1桁の単位=128ビット
前回の例６（P=2189*3 + 1、1桁の単位=64ビット)は192ビット乗算を行いますが、例８では160ビットで済むのでその分若干早くなりました。
同様に例７（P=2316*13 + 1、1桁の単位=128ビット)は320ビット乗算ですが、例９では288ビットで済みます。
例８、９で、64ビット浮動小数点計算の複素ＦＦＴと大体同じ位の速度になりました。
現状ではここ辺りが限界でしょうか。
円周率の計算時間は、サラミンとブレントの改良公式はチュドノフスキーの公式での２倍強かかります。
チュドノフスキーの公式の計算では、２５０万桁が約５０秒、．．．、４０００万桁が約２５分、８０００万桁が１時間弱かかりました。これにはＫＡＲＡＴＳＵＢＡ方式も併用（小さい桁の乗算で）しています。
チュドノフスキーの公式では殆どの計算が乗算なので、乗算速度が上がらないと性能が向上しません。しかし算術幾何平均を使った公式群では冪乗根の計算が主体なので、乗算速度が同じでもニュートン法の計算手順のなかでの改善で、見直せる事があるかもしれません（でも２倍の差は埋まらないでしょう、多分）。
しかしその前に計算結果の出力（２進１０進変換）処理の高速化が必要になりました。こちらは全然性能など考えていませんでしたが、桁数が上がるにつれて２進１０進変換時間の方が長くかかる様になってしまいました。なので前記桁数は計算確定の時間ではありませんので、悪しからず。
サラミンとブレントの改良公式で１０００万桁、チュドノフスキーの公式で５００万桁までは、なんとか出力して確認できました。
他に、ワラを掴む的な準備（下調べ）として、Toom-Cook法（４分割）の計算を確認してみました。参考記事では行列表記でさらりと説明されていますが、行列（式）は得意ではないので普通の式に置き換えて、十進ＢＡＳＩＣで計算してみました。すると計算結果が異なる場合があり、正誤のいくつかを見比べて、どうも係数が違っている（らしい？）事が分かりました。正解はｘｘ元連立方程式を解かないと出てこないらしいので、それは難しそうだからスキップして、正しい計算結果となるような係数を逆算で推定してみました（”７２”→”１８００”）。
性能は度外視して、理論通り？と思われる式で計算した結果が次のとおりです。KARATSUBA方式の倍に分割しただけなのに、なんでこんなに式が多くなるんだか．．．。
まだやれそうな（やらなきゃならない）事があるので、取敢えず良し々々。
≪プラグラム例と結果≫
（１６桁×１６桁を４桁ずつに分解して計算しています）
OPTION ARITHMETIC RATIONAL
LET dn=4
LET d0=10^4
LET xx=9988776655443322
LET yy=9999222288881111
PRINT "xx =";xx,"yy =";yy
PRINT "非分割乗算 xx * yy =";
PRINT USING REPEAT$("-",8*dn-1) & "%":xx*yy
LET x0=MOD(xx,d0)
LET x1=MOD(INT(xx/d0),d0)
LET x2=MOD(INT(xx/d0^2),d0)
LET x3=MOD(INT(xx/d0^3),d0)
LET y0=MOD(yy,d0)
LET y1=MOD(INT(yy/d0),d0)
LET y2=MOD(INT(yy/d0^2),d0)
LET y3=MOD(INT(yy/d0^3),d0)
PRINT "xx=";x3;x2;x1;x0,"yy=";y3;y2;y1;y0
LET u6= -x3 +2*x2 -4*x1 +8*x0
LET u5=-8*x3 +4*x2 -2*x1 +x0
LET u4= -x3 +x2 -x1 +x0
LET u3= x0
LET u2= x3 +x2 +x1 +x0
LET u1= 8*x3 +4*x2 +2*x1 +x0
LET u0= x3
LET v6= -y3 +2*y2 -4*y1 +8*y0
LET v5=-8*y3 +4*y2 -2*y1 +y0
LET v4= -y3 +y2 -y1 +y0
LET v3= y0
LET v2= y3 +y2 +y1 +y0
LET v1= 8*y3 +4*y2 +2*y1 +y0
LET v0= y3
PRINT "u=";u6;u5;u4;u3;u2;u1;u0
PRINT "v=";v6;v5;v4;v3;v2;v1;v0
LET uv6=u6*v6
LET uv5=u5*v5
LET uv4=u4*v4
LET uv3=u3*v3
LET uv2=u2*v2
LET uv1=u1*v1
LET uv0=u0*v0
PRINT "uv=";uv6;uv5;uv4;uv3;uv2;uv1;uv0
LET w6= 360*uv0
LET w5= -4*uv6-10*uv5+120*uv4+180*uv3 -40*uv2 +6*uv1 +180*uv0
LET w4= 15*uv5 -60*uv4 +90*uv3 -60*uv2 +15*uv1-1800*uv0
LET w3= 20*uv6+20*uv5-540*uv4-900*uv3+140*uv2 -900*uv0
LET w2= -15*uv5+240*uv4-450*uv3+240*uv2 -15*uv1+1440*uv0
LET w1=-16*uv6-10*uv5+240*uv4+720*uv3 +80*uv2 -6*uv1 +720*uv0
LET w0= 360*uv3
LET w6=w6/360
LET w5=w5/360
LET w4=w4/360
LET w3=w3/360
LET w2=w2/360
LET w1=w1/360
LET w0=w0/360
PRINT "w=";w6;w5;w4;w3;w2;w1;w0
LET w1=w1+INT(w0/d0)
LET w0=MOD(w0,d0)
LET w2=w2+INT(w1/d0)
LET w1=MOD(w1,d0)
LET w3=w3+INT(w2/d0)
LET w2=MOD(w2,d0)
LET w4=w4+INT(w3/d0)
LET w3=MOD(w3,d0)
LET w5=w5+INT(w4/d0)
LET w4=MOD(w4,d0)
LET w6=w6+INT(w5/d0)
LET w5=MOD(w5,d0)
LET w7=w7+INT(w6/d0)
LET w6=MOD(w6,d0)
LET fmt1$=REPEAT$("-",dn-1) & "%"
LET fmt2$=REPEAT$("%",dn)
PRINT "Toom-Cook法(４分割) zz = ";
PRINT USING fmt1$:w7;
PRINT USING fmt2$:w6;
PRINT USING fmt2$:w5;
PRINT USING fmt2$:w4;
PRINT USING fmt2$:w3;
PRINT USING fmt2$:w2;
PRINT USING fmt2$:w1;
PRINT USING fmt2$:w0;
END
xx = 9988776655443322 yy = 9999222288881111
非分割乗算 xx * yy =99879998171764182850815056890742
xx= 9988 7766 5544 3322 yy= 9999 2222 8888 1111
u= 9944 -56606 -4444 3322 26620 125378 9988
v=-32219 -87769 -15554 1111 22220 107767 9999
uv=-320385736 4968252014 69121976 3690742 591496400 13511610926 99870012
w= 99870012 99845570 161463852 125656322 65284582 35685320 3690742
Toom-Cook法(４分割) zz = 99879998171764182850815056890742
（注）本文は、数学素人の個人的な記録（単なるメモ）なので、誤りがある事などをご了解ください。
≪参考文献・記事≫
計算代数（加古富志雄H26.4.24）
円周率．ｊｐ （ウィキペディア）