円周率の計算（続き： ＦＦＴなど乗算方法）
（2015.4.26作成）
桁数の大きい数値の演算（乗算など）は、CPUの1命令や普通の言語での1ステートメントで実行することはできません（十進BASICの有理数モードでは、相当の桁数まで同一のプログラムで計算できます）。
最も単純な方法は、個々の変数に2～10桁位ずつを入れて筆算で行うような手順をプログラムすることです。しかし桁数が大きくなっていくと計算量も膨大になっていきます。100桁の乗算を筆算で行う場合、九九を10000回唱えなければなりません。桁数が10倍になると計算回数（乗算）が100倍になります。
（１）KARATSUBAの方式
それで頭の良い人、KARATSUBA（カラツバ）と言う人、が考えてくれました。2桁同士の掛け算で普通なら九九が4回必要になります。
（10*A + B）*(10*C + D) ⇒ 100*A*C + 10*（A*D + B*C) + B*D
しかし一寸捻って、
（10*A + B）*(10*C + D) ⇒ 100*A*C + 10*（(A - B)*(D - C) + A*C + B*D) + B*D
こんな具合に、3回の九九（掛け算）だけで済んでしまいます。ただし加算が1回、減算が2回増えてしまいますが。
人間が筆算をする場合には、こんな入り組んだ計算をする方が却って手間取りますが、計算機で計算をする場合は違ってきます（プログラム作成が手間取るだけです）。計算機では、乗算は加減算の10～20倍くらい時間がかかりますので、乗算回数の量が全体の処理時間に対する大きな比重を占めます。乗算回数の比率が全体の計算時間の比率を左右するのです（当然、ディスクの読み書きなどは除外しての話です）。
桁数が大きい場合、半分の桁数になった掛け算（上記のA＊C、B＊Dなど）に同じ手順を適用します。100桁の場合、最後（1桁）まで半分々々に出来ませんが、理論的（計算上）には、1500回未満の掛け算（九九）で済む事になります。
とにかく考え方は誰でも直ぐに理解できますし実現も簡単なので、試す価値があります。
より高性能な乗算が出来るFFT（後述）では大体が浮動小数点計算を利用をします。計算機の浮動小数点計算と整数計算の処理時間（速度）を比べると、PCに使われてるCPUの種類によって差があります。
浮動小数点計算速度 ≪ 整数計算速度とか、浮動小数点計算速度 ＜ 整数計算速度という違いです。またFFTを乗算で使う場合、付随する計算がかなりあって、計算桁数がある程度大きくならないと効果が表れてきません。上の様に浮動小数点と整数の計算速度も絡んできます。1兆桁などでは勝負になりませんが、1000～10000桁位までは案外KARATSUBA方式の方が性能が高くなります。
KARATSUBA方式は2分割で加減算を工夫して乗算を削減していますが、同様の考え方で3分割、4分割、．．．などで乗算回数を減らす方式があるようです。しかし分割数が増えると変形の式がどんどん複雑怪奇？になって、なかなか試そうという意欲が湧いてこなくなります。
（２）FFT
こちらはKARATSUBA方式のように一言で説明できませんし、細かく（理論的に）語る知識も能力もありません。多数桁乗算の関連だけでザックリの話をすると、次の様になります。
まずはフーリエ変換、これは数学的研究や物理現象の分析などに利用されるようですが、関数にある種の操作を加えて別の関数に変換（⇒）するわけです。そして元の関数に戻す変換（➡）もあります（逆変換）。
y＝f（x） ⇒ y＝F（x）、 y＝F（x） ➡ y＝f（x）
関数なので、x、yは概ね連続値ですが、計算機でデータ（音の信号など）処理を行う場合、一定時間ごとの独立した数値などを扱うことになります。それで離散フーリエ変換（と逆変換）が登場します。離散（数値列）なので変換手順は異なりますが、考え方は一緒です。（連続関数での積分 → 数値列同士の積の総和（シグマ）計算、という感じ）。xを等間隔であるとして省略すると、
(y0, y1, . . .) ⇒ （Y0, Y1, . . .)、 （Y0, Y1, . . .) ➡ (y0, y1, . . .)
このフーリエ変換（連続、離散とも）、畳込演算という便利な操作が出来ます。ここからは離散（数列）だけで話をすすめます。今2桁数値の掛け算を考えます。
(10*x1 + x2)*(10*y1 + y2)
筆算の場合、
x1*y1=10*p1 + p2、 x1*y2=10*q1 + q2、
x2*y1=10*r1 + r2、 x2*y2=10*s1 + s2
A=100*(10*p1 + p2) + 10*(10*q1 + q2 + 10*r1 + r2) + 10*s1 + s2
=1000*p1 + 100*(p2 + q1 + r1) + 10*(q2 + r2 + s2) + s2
の4桁数値になります。
これをフーリエ変換（結果が4桁になるので、通しで4桁変換を行います）、畳込演算、フーリエ逆変換で計算（最初のｘ、ｙは虚数をゼロとして、全てを複素数計算）すると、
フーリエ変換： (0、0、x1、x2)⇒(u1、u2、u3、u4)、
(0、0、y1、y2)⇒(v1、v2、v3、v4)
畳込演算： w1=u1*v1、 w2=u2*v2、 w3=u3*v3、 w4=u4*v4
フーリエ逆変換： (w1、w2、w3、w4)➡(z1、z2、z3、z4)
// z1=0、 z2=x1*y1、 z3=x1*y2 + x2*y1、 z4=x2*y2 //
後処理： z4からz1まで、筆算のように桁上げ計算をする必要があります（ここ
桁数の大きい数値の演算（乗算など）は、CPUの1命令や普通の言語での1ステートメントで実行することはできません（十進BASICの有理数モードでは、相当の桁数まで同一のプログラムで計算できます）。
最も単純な方法は、個々の変数に2～10桁位ずつを入れて筆算で行うような手順をプログラムすることです。しかし桁数が大きくなっていくと計算量も膨大になっていきます。100桁の乗算を筆算で行う場合、九九を10000回唱えなければなりません。桁数が10倍になると計算回数（乗算）が100倍になります。
（１）KARATSUBAの方式
それで頭の良い人、KARATSUBA（カラツバ）と言う人、が考えてくれました。2桁同士の掛け算で普通なら九九が4回必要になります。
（10*A + B）*(10*C + D) ⇒ 100*A*C + 10*（A*D + B*C) + B*D
しかし一寸捻って、
（10*A + B）*(10*C + D) ⇒ 100*A*C + 10*（(A - B)*(D - C) + A*C + B*D) + B*D
こんな具合に、3回の九九（掛け算）だけで済んでしまいます。ただし加算が1回、減算が2回増えてしまいますが。
人間が筆算をする場合には、こんな入り組んだ計算をする方が却って手間取りますが、計算機で計算をする場合は違ってきます（プログラム作成が手間取るだけです）。計算機では、乗算は加減算の10～20倍くらい時間がかかりますので、乗算回数の量が全体の処理時間に対する大きな比重を占めます。乗算回数の比率が全体の計算時間の比率を左右するのです（当然、ディスクの読み書きなどは除外しての話です）。
桁数が大きい場合、半分の桁数になった掛け算（上記のA＊C、B＊Dなど）に同じ手順を適用します。100桁の場合、最後（1桁）まで半分々々に出来ませんが、理論的（計算上）には、1500回未満の掛け算（九九）で済む事になります。
とにかく考え方は誰でも直ぐに理解できますし実現も簡単なので、試す価値があります。
より高性能な乗算が出来るFFT（後述）では大体が浮動小数点計算を利用をします。計算機の浮動小数点計算と整数計算の処理時間（速度）を比べると、PCに使われてるCPUの種類によって差があります。
浮動小数点計算速度 ≪ 整数計算速度とか、浮動小数点計算速度 ＜ 整数計算速度という違いです。またFFTを乗算で使う場合、付随する計算がかなりあって、計算桁数がある程度大きくならないと効果が表れてきません。上の様に浮動小数点と整数の計算速度も絡んできます。1兆桁などでは勝負になりませんが、1000～10000桁位までは案外KARATSUBA方式の方が性能が高くなります。
KARATSUBA方式は2分割で加減算を工夫して乗算を削減していますが、同様の考え方で3分割、4分割、．．．などで乗算回数を減らす方式があるようです。しかし分割数が増えると変形の式がどんどん複雑怪奇？になって、なかなか試そうという意欲が湧いてこなくなります。
（２）FFT
こちらはKARATSUBA方式のように一言で説明できませんし、細かく（理論的に）語る知識も能力もありません。多数桁乗算の関連だけでザックリの話をすると、次の様になります。
まずはフーリエ変換、これは数学的研究や物理現象の分析などに利用されるようですが、関数にある種の操作を加えて別の関数に変換（⇒）するわけです。そして元の関数に戻す変換（➡）もあります（逆変換）。
y＝f（x） ⇒ y＝F（x）、 y＝F（x） ➡ y＝f（x）
関数なので、x、yは概ね連続値ですが、計算機でデータ（音の信号など）処理を行う場合、一定時間ごとの独立した数値などを扱うことになります。それで離散フーリエ変換（と逆変換）が登場します。離散（数値列）なので変換手順は異なりますが、考え方は一緒です。（連続関数での積分 → 数値列同士の積の総和（シグマ）計算、という感じ）。xを等間隔であるとして省略すると、
(y0, y1, . . .) ⇒ （Y0, Y1, . . .)、 （Y0, Y1, . . .) ➡ (y0, y1, . . .)
このフーリエ変換（連続、離散とも）、畳込演算という便利な操作が出来ます。ここからは離散（数列）だけで話をすすめます。今2桁数値の掛け算を考えます。
(10*x1 + x2)*(10*y1 + y2)
筆算の場合、
x1*y1=10*p1 + p2、 x1*y2=10*q1 + q2、
x2*y1=10*r1 + r2、 x2*y2=10*s1 + s2
A=100*(10*p1 + p2) + 10*(10*q1 + q2 + 10*r1 + r2) + 10*s1 + s2
=1000*p1 + 100*(p2 + q1 + r1) + 10*(q2 + r2 + s2) + s2
の4桁数値になります。
これをフーリエ変換（結果が4桁になるので、通しで4桁変換を行います）、畳込演算、フーリエ逆変換で計算（最初のｘ、ｙは虚数をゼロとして、全てを複素数計算）すると、
フーリエ変換： (0、0、x1、x2)⇒(u1、u2、u3、u4)、
(0、0、y1、y2)⇒(v1、v2、v3、v4)
畳込演算： w1=u1*v1、 w2=u2*v2、 w3=u3*v3、 w4=u4*v4
フーリエ逆変換： (w1、w2、w3、w4)➡(z1、z2、z3、z4)
// z1=0、 z2=x1*y1、 z3=x1*y2 + x2*y1、 z4=x2*y2 //
後処理： z4からz1まで、筆算のように桁上げ計算をする必要があります（ここ
での浮動小数点計算は整数部分だけを使います）
例えば、各変数が2桁処理でx1=85、y1=31の場合、z2=2635、それを
z1=26、z2=35に桁上げ処理するという具合です
となるのです。これで掛け算が出来た事になります。2桁では掛け算回数（4回）が減っていませんが、100桁の場合を考えると、筆算では10000回、フーリエ変換と畳込演算を使うと100回で済むという驚異の差になります。
問題はフーリエ変換と逆変換です。ここでも実は掛け算が必要なのですが、100桁の場合、10000回ずつになります。なのでフーリエ変換をそのまま使うと全然高速化にはなりません。
そしてまた頭の良い人がFFT（高速フーリエ変換）を考えだしてくれました。
フーリエ変換と逆変換は、x1やx2にsinやcosの値を掛けて合計するという処理になっています。4個の数列の場合、
X1=sin(0)*x1 + sin(0)*x2 + sin(0)*x3 + sin(0)*x4
X2=sin(0)*x1 + sin(π/2)*x2 + sin(π)*x3 + sin(3π/2)*x4
X3=sin(0)*x1 + sin(π)*x2 + sin(2π)*x3 + sin(3π)*x4
X4=sin(0)*x1 + sin(3π/2)*x2 + sin(3π)*x3 + sin(9π/2)*x4
大体こんな感じです（雰囲気を掴んでもらうための式なので、正確ではありません）。4桁の場合掛け算が16回必要です。しかしπの有理数倍のsinやcosは、お互いに同値だったり符号が違うだけだったり何らかの関係があって、16回丸々計算しなくても済むと気づいた訳です。結果、
x →（FFT）→ X、 （”FFT” ≡ ”⇒”）
なるフーリエ変換と同等で、必要な掛け算回数が桁数に比例する程度で済む変換方法が編出されました。これが世に言う高速フーリエ変換（FFT）です。これで初めてFFTと畳込演算を組合わせて、筆算よりはるかに高速な掛け算が可能となりました。
必要な計算精度： 上の例でフーリエ逆変換で出力されたz3の内訳（式: 桁上げ処理前）を見ると、z3の値は最大で
問題はフーリエ変換と逆変換です。ここでも実は掛け算が必要なのですが、100桁の場合、10000回ずつになります。なのでフーリエ変換をそのまま使うと全然高速化にはなりません。
そしてまた頭の良い人がFFT（高速フーリエ変換）を考えだしてくれました。
フーリエ変換と逆変換は、x1やx2にsinやcosの値を掛けて合計するという処理になっています。4個の数列の場合、
X1=sin(0)*x1 + sin(0)*x2 + sin(0)*x3 + sin(0)*x4
X2=sin(0)*x1 + sin(π/2)*x2 + sin(π)*x3 + sin(3π/2)*x4
X3=sin(0)*x1 + sin(π)*x2 + sin(2π)*x3 + sin(3π)*x4
X4=sin(0)*x1 + sin(3π/2)*x2 + sin(3π)*x3 + sin(9π/2)*x4
大体こんな感じです（雰囲気を掴んでもらうための式なので、正確ではありません）。4桁の場合掛け算が16回必要です。しかしπの有理数倍のsinやcosは、お互いに同値だったり符号が違うだけだったり何らかの関係があって、16回丸々計算しなくても済むと気づいた訳です。結果、
x →（FFT）→ X、 （”FFT” ≡ ”⇒”）
なるフーリエ変換と同等で、必要な掛け算回数が桁数に比例する程度で済む変換方法が編出されました。これが世に言う高速フーリエ変換（FFT）です。これで初めてFFTと畳込演算を組合わせて、筆算よりはるかに高速な掛け算が可能となりました。
必要な計算精度： 上の例でフーリエ逆変換で出力されたz3の内訳（式: 桁上げ処理前）を見ると、z3の値は最大で
最大桁数 ＊ （1桁の最大数） ＊ （1桁の最大数）
の値を桁落ちなく保持できる必要がある事が分かります。PCのCPUで倍精度浮動小数点は52ビット（機構的には64ビット）の精度がありますので、1桁を2進16ビット（0～65535）とすれば、理論上処理可能な最大桁数は、約100万桁になります。実際は途中計算で誤差が入ってくるので、1ビット分位は目減りする（1／2になる）可能性が高いです（つまり50万桁程）。2進16ビットは10進4.8桁分になるので、10進で250万桁くらいの掛け算ができる訳です。同様に1桁を2進8ビットとすれば、普通のPCではほぼ無限大まで計算できそうです。ただFFT処理で1桁毎かなりのオーバヘッドがかかります。1桁を16ビットで処理しても8ビットで処理しても大体同じくらいのオーバヘッドです。つまりは8ビット処理にすると、処理速度が半減してしまうという事です。
オーバヘッドは速度だけでなく、メモリにも絡みます。通常のフーリエ変換（FFTを含む）は、入力、内部データ、出力ともに複素数になります。入力時は虚数には0を入れ、出力時は虚数（0の筈ですが）を無視します。8ビット処理でも、FFT処理としては倍精度浮動小数点データ（64ビット）が2個使われることになります。
性能改善：
（ａ）ハード命令の精度をフル活用
直ぐに思いついたのは、CPUの機構として持っている仮数64ビットのデータ処理です。それを利用すれば、1桁を16ビットとしても理論上は40億桁まで処理できる勘定です。しかし手持ちのc言語ではサポートしていません（昔の16ビット版ではサポートしていたそうですが）。自前のアセンブルコードで書かなければならずずっと放ったままでおきましたが、最近漸くヤル気を出して完成しました。8ビット処理時の倍の性能になるかと期待したのですが、さにあらず。いくらか早い程度にしかならず、少々ガッカリしました。自前のアセンブラコードがコンパイラのよりかなり劣るということなのかな？メモリアクセスの単位も変わってくる（64ビット（52ビット仮数＋指数）と80ビット（64ビット仮数＋指数）の違い）ので一概にコードだけの差かどうかは不明です。追及はまたいつかという事にしておきます。
（ｂ）入出力が実数（実際は整数だけ）に限られるという条件の活用
FFTの使用目的が多桁の乗算なので、入出力で虚数は使っていません（不要）。それで普通なら乗数と被乗数を別々にFFTするところを、乗数を実数、被乗数を虚数とする複素数を入力として1回のFFTを行い、変換後に2つの変数に分離してから畳込演算を行います。逆変換では結果の1変数を処理するわけですが、その場合事前にうまい演算を施す事により、半分の桁数の逆変換で、答えとなる数値列が実数部と虚数部交互に出力されるようにできます。これらの対策で処理速度が何％か向上します。
イメージで書くと次の様になります。通常の場合
FFT： (0、0、x1、x2)⇒(u1、u2、u3、u4)、
(0、0、y1、y2)⇒(v1、v2、v3、v4)
畳込演算： w1=u1*v1、 w2=u2*v2、 w3=u3*v3、 w4=u4*v4
逆FFT： (w1、w2、w3、w4)➡(z1、z2、z3、z4)
z1=0、 z2=x1*y1、 z3=x1*y2 + x2*y1、 z4=x2*y2
改善対策の場合
FFT： (0、0、x1+iy1、x2+iy2)⇒(uv1、uv2、uv3、uv4)
分離計算： (uv1、uv2、uv3、uv4) → (u1、u2、u3、u4)、 (v1、v2、v3、v4)
畳込演算： w1=u1*v1、 w2=u2*v2、 w3=u3*v3、 w4=u4*v4
圧縮計算： (w1、w2、w3、w4) → (-、-、t1、t2)
逆FFT： (t1、t2)➡(z1+iz2、z3+iz4）
// z1=0、 z2=x1*y1、 z3=x1*y2 + x2*y1、 z4=x2*y2 //
分離・圧縮計算はFFT計算より更に訳分からないまま見様見真似で作成しましたが、何とか作れました。
（ｃ）浮動小数点ビット数の拡大
浮動小数点計算、処理単位を増やしてオーバヘッドの比率を下げられれば、性能向上するのではないかと思い、
ソフトでハード命令64ビット精度を超える浮動小数点計算を作る事にしました。最初は仮数128ビットで処理を32ビット／桁に。しかしソフトで動かすとかなり遅い。10倍まではいかないにしても、全然改善になりません。それで悪あがきに仮数192ビットの64ビット／桁処理版を作ってみました。前作よりましになりましたが、ハード命令を使うのにはかなわず、この方向の改善？は止めにしました。
（ｄ）整数環FFT
フーリエ変換は、もともとが複素数の領域で定義されていますが、高速乗算に使うだけなら、フーリエ変換本来の効果（意味）を無視して、整数領域だけで行えるという説明があったのを見て、それを試す事にしました。
或る値（P）を法とする整数環のなかでFFTと同等の演算が行なえる（らしい）と。それをFFTと呼べるのか私の知見でコメント出来ませんし、私にとっては説明が簡単（簡明）過ぎて理解しきれず、手探りで進むしかありませんでした。
計算は整数に限られるので、通常（複素FFT）の計算式から虚数の項を取り除いて、変換で使う式を推定します。加減乗算は全てモジュラー計算とします。
或る値Pは、かなりデカイ数値（>264 ）が必要となりそうなでのモジュラー計算（剰余計算）をどうしようか? ネット検索して、モンゴメリー乗算（乗算のみで、剰余計算（モジュラー計算）をしてしまうという技）なる計算手法を見つけました。
肝心なPはどうするか? 任意の値で整数環になるとは思えないし（否、整数環となってもFFTとして使えないか、…）、ここは素人考えで素数なら望ましい整数環になるんじゃなかろーか．．．などと勝手に決めつけました。
ではデカイ素数（値）をどうやってみつけるのか？ネットでミラー・ラビン素数判定法というのを見つけたので、十進BASICで試してみました。試しながらも、つらつらと考えました。この方法は被検査数値が素数でしたという確定結果を答えてくれません（たぶん素数でしょう、．．．みたいな。＝＝＞m(_ _)mすみません、決してこのアルゴリズムの批判とかでなく、素人の素朴な感想です。m(_ _)m）。それで他に簡単で確定結果を出してくれる方法は無いのか? 仕方なく、書店で探すことにしました。なるべく本の在庫が豊富そうな、池袋ジュンク堂書店へ。有りました !! 「[講座]数学の考え方⑯初等整数論」木田祐司（著）、この中に＜原始根テスト＞なる方法が書かれていました。
＜原始根テスト＞
或る数nが素数か？
＝＞n - 1を素因数分解する。 （例えば n - 1 = a1k1 * a2k2 * a3k3 だったとします）
＝＞その場合に或る数qを仮定して、nと各因数(a1、a2、a3)とqで作られた判定式をクリアしたら、
＝＞qはnの原始根（と言うらしい）で、nは素数である。
という感じです。ｑをどこまで探せば良いのか分かりませんが、時間がかかるので、2～1000に絞って、それで駄目なら検査終了（素数ではなさそう）とすることにしました。n －1の素因数分解は巨大な数では難しい訳ですが、別に暗号に使う為ではないので、n－1の因数が、判りやすく且つ種類が多くならないような数を候補として選びます。
これも十進BASICでササッと作成。探すのは、2(32＊i)より小さくて、出来るだけそれに近い素数。さらに都合の良い1の累乗根(回転子：複素数版のsin／cosに相当）が見つかるように、
2((32*i) - j) * k + 1、 2 <= k < 2j
の形の数値で上記原始根テストを行います。
素数が見つかったら回転子（整数環での1の累乗根）を求める訳ですが、数値がデカイですし何より累乗根を計算できないので、逆に数値を適当に選んでは累乗してみて1になるまでの限界を見極めて選びます。（全数の累乗はとても計算できないので、22、23、．．．、2xの各数値で累乗して、1になる限界をチェックします）。前述の様な素数が見つかった場合、大体2(32*i - j)乗で1になるω（回転子？）があるみたいです。
P=2u * v + 1
1≡ωw mod(P) (w=2u または w=2(u - 1))
というP、u、v、ωを探した訳です。その結果見つかった手頃なもので、整数環FFTを作って性能比較をしてみました。
64～320ビット同士の掛け算やモンゴメリ乗算は、大体c言語のインラインアセンブルコードを主体として作りました。
例１．P=228*13 + 1 (16進表示: 0xD0000001) (10進表示: 3489660929)
ω=11 (w=228)
1桁の単位=8ビット
例２．P=230*3 + 1 (16進表示: 0xC0000001) (10進表示: 3221225473)
ω=13 (w=230)
1桁の単位=8ビット
例３．P=259*27 + 1 (16進表示: 0xD80000000, 1)
(10進表示: 15564440312192434177)
ω=87 (w=259)
1桁の単位=16ビット
例４．P=292*7 + 1 (16進表示: 0x70000000, 0, 1)
(10進表示: 34662321099990647697175478273)
ω=11 (w=292)
1桁の単位=32ビット
例５．P=2121*81 + 1 (16進表示: 0xA2000000, 0, 0, 1)
(10進表示: 215334935317156371410416743765415821313)
ω=57 (w=2120)
1桁の単位=32ビット
例６．P=2189*3 + 1 (16進表示: 0x60000000, 0, 0, 0, 0, 1)
(10進表示: 2353913150770005286438421033702874906038383291674012942337)
ω=3 (w=2188)
1桁の単位=64ビット
例７．P=2316*13 + 1 (16進表示: 0xD0000000, 0, 0, 0, 0, 0, 0, 0, 0, 1)
(10進表示: 1735489466685739441945955136262761093114697424414780375581971306355553527196770446893656695635969)
ω=7 (w=2316)
1桁の単位=128ビット
性能比較：
Pが小さいと桁当たりの処理オーバヘッドが大きくなり、Pが大きいと手作り掛け算処理の速度が影響してきます。性能（と言ってもどんぐりの背比べです）的には、例４）P＝292*7+1、例６）P＝2189*3+1、例７）P＝2316 *13+1 がまあまあでした。浮動小数点（複素数）FFTといい勝負です。
KARATSUBAのところで書いたように、CPU（整数と浮動小数点の演算速度の違い）によってこの差も異なります。
今使ってるノートPCは、
（PC1）HP1000-1401（Win8、CPUはIntel Celeron CPU B830 CLK1.8GHz）
で、その前にメインで使っていたのは、
（PC2）ディスクトップＰＣ（ＰＣショップオリジナル）
（WinXP、CPUはAMD Athlon 64X2 Dual Core Proc 5400+ CLK2.8GHz）
ですが、PC2はPC1に比べて浮動小数点計算が相対的に遅いです。
いろいろ作った乗算処理（FFTとKARATSUBA方式）の速度（経過時間：秒）は、現在値としてこのくらいです。
我が家のPCでは、メモリ食いの浮動小数点（複素数）FFTを使って20million桁以上の計算は実行できません。
メモリ効率をみると、整数環FFTのほうが大勝?です。8ビットを処理するためにFFT処理で何ビット必要かを比べると、
標準（c言語の倍精度浮動小数点52ビット精度）： 128ビット（十進300万桁以上の計算の場合）
64ビット精度浮動小数点： 96ビット
手作り192ビット精度浮動小数点： 48ビット
整数環FFT（例４）： 24ビット
整数環FFT（例６）： 24ビット
整数環FFT（例７）： 20ビット
これから：
なにか改善?の気力が湧いて実行できたら、また記録を残したいと思います。
追記：
整数環FFT（一般のも含めて?）の回転子、入力桁の数を分離できるだけの値（桁数の累乗根）で良いのでしょうか?ヤル気が出たら調べてみましょう。しかし計算はちゃんとできてるようなので、今のところは問題ありません。この辺り、素人なので無駄なところに労力（神経）を割いてしまっているような気がします。
前の記事で、ラマヌジャン公式の計算ミスについて書きましたが、調べたら大事ではなく解決しました。最終式では”B”値を使っていない（チュドノフスキーの公式も同様）ので、最後の1項になるときの計算で、時間節約のためにB計算をしないようにしていたのですが、その判定条件を誤っていただけでした。よかった、よかった。
（注）本文は、数学素人の個人的な記録（単なるメモ）なので、誤りがある事などをご了解ください。
≪参考文献・記事≫
計算機代数入門（野呂正行）
数値計算の基礎（川上一郎）
「[講座]数学の考え方（16）初等整数論」木田祐司（著）朝倉書店
円周率.jp （ウィキペディア）
モンゴメリー乗算 （ウィキペディア）
FFT（高速フーリエ・コサイン・サイン変換）の概略と設計法
（Takuya OOURA）
の値を桁落ちなく保持できる必要がある事が分かります。PCのCPUで倍精度浮動小数点は52ビット（機構的には64ビット）の精度がありますので、1桁を2進16ビット（0～65535）とすれば、理論上処理可能な最大桁数は、約100万桁になります。実際は途中計算で誤差が入ってくるので、1ビット分位は目減りする（1／2になる）可能性が高いです（つまり50万桁程）。2進16ビットは10進4.8桁分になるので、10進で250万桁くらいの掛け算ができる訳です。同様に1桁を2進8ビットとすれば、普通のPCではほぼ無限大まで計算できそうです。ただFFT処理で1桁毎かなりのオーバヘッドがかかります。1桁を16ビットで処理しても8ビットで処理しても大体同じくらいのオーバヘッドです。つまりは8ビット処理にすると、処理速度が半減してしまうという事です。
オーバヘッドは速度だけでなく、メモリにも絡みます。通常のフーリエ変換（FFTを含む）は、入力、内部データ、出力ともに複素数になります。入力時は虚数には0を入れ、出力時は虚数（0の筈ですが）を無視します。8ビット処理でも、FFT処理としては倍精度浮動小数点データ（64ビット）が2個使われることになります。
性能改善：
（ａ）ハード命令の精度をフル活用
直ぐに思いついたのは、CPUの機構として持っている仮数64ビットのデータ処理です。それを利用すれば、1桁を16ビットとしても理論上は40億桁まで処理できる勘定です。しかし手持ちのc言語ではサポートしていません（昔の16ビット版ではサポートしていたそうですが）。自前のアセンブルコードで書かなければならずずっと放ったままでおきましたが、最近漸くヤル気を出して完成しました。8ビット処理時の倍の性能になるかと期待したのですが、さにあらず。いくらか早い程度にしかならず、少々ガッカリしました。自前のアセンブラコードがコンパイラのよりかなり劣るということなのかな？メモリアクセスの単位も変わってくる（64ビット（52ビット仮数＋指数）と80ビット（64ビット仮数＋指数）の違い）ので一概にコードだけの差かどうかは不明です。追及はまたいつかという事にしておきます。
（ｂ）入出力が実数（実際は整数だけ）に限られるという条件の活用
FFTの使用目的が多桁の乗算なので、入出力で虚数は使っていません（不要）。それで普通なら乗数と被乗数を別々にFFTするところを、乗数を実数、被乗数を虚数とする複素数を入力として1回のFFTを行い、変換後に2つの変数に分離してから畳込演算を行います。逆変換では結果の1変数を処理するわけですが、その場合事前にうまい演算を施す事により、半分の桁数の逆変換で、答えとなる数値列が実数部と虚数部交互に出力されるようにできます。これらの対策で処理速度が何％か向上します。
イメージで書くと次の様になります。通常の場合
FFT： (0、0、x1、x2)⇒(u1、u2、u3、u4)、
(0、0、y1、y2)⇒(v1、v2、v3、v4)
畳込演算： w1=u1*v1、 w2=u2*v2、 w3=u3*v3、 w4=u4*v4
逆FFT： (w1、w2、w3、w4)➡(z1、z2、z3、z4)
z1=0、 z2=x1*y1、 z3=x1*y2 + x2*y1、 z4=x2*y2
改善対策の場合
FFT： (0、0、x1+iy1、x2+iy2)⇒(uv1、uv2、uv3、uv4)
分離計算： (uv1、uv2、uv3、uv4) → (u1、u2、u3、u4)、 (v1、v2、v3、v4)
畳込演算： w1=u1*v1、 w2=u2*v2、 w3=u3*v3、 w4=u4*v4
圧縮計算： (w1、w2、w3、w4) → (-、-、t1、t2)
逆FFT： (t1、t2)➡(z1+iz2、z3+iz4）
// z1=0、 z2=x1*y1、 z3=x1*y2 + x2*y1、 z4=x2*y2 //
分離・圧縮計算はFFT計算より更に訳分からないまま見様見真似で作成しましたが、何とか作れました。
（ｃ）浮動小数点ビット数の拡大
浮動小数点計算、処理単位を増やしてオーバヘッドの比率を下げられれば、性能向上するのではないかと思い、
ソフトでハード命令64ビット精度を超える浮動小数点計算を作る事にしました。最初は仮数128ビットで処理を32ビット／桁に。しかしソフトで動かすとかなり遅い。10倍まではいかないにしても、全然改善になりません。それで悪あがきに仮数192ビットの64ビット／桁処理版を作ってみました。前作よりましになりましたが、ハード命令を使うのにはかなわず、この方向の改善？は止めにしました。
（ｄ）整数環FFT
フーリエ変換は、もともとが複素数の領域で定義されていますが、高速乗算に使うだけなら、フーリエ変換本来の効果（意味）を無視して、整数領域だけで行えるという説明があったのを見て、それを試す事にしました。
或る値（P）を法とする整数環のなかでFFTと同等の演算が行なえる（らしい）と。それをFFTと呼べるのか私の知見でコメント出来ませんし、私にとっては説明が簡単（簡明）過ぎて理解しきれず、手探りで進むしかありませんでした。
計算は整数に限られるので、通常（複素FFT）の計算式から虚数の項を取り除いて、変換で使う式を推定します。加減乗算は全てモジュラー計算とします。
或る値Pは、かなりデカイ数値（>264 ）が必要となりそうなでのモジュラー計算（剰余計算）をどうしようか? ネット検索して、モンゴメリー乗算（乗算のみで、剰余計算（モジュラー計算）をしてしまうという技）なる計算手法を見つけました。
肝心なPはどうするか? 任意の値で整数環になるとは思えないし（否、整数環となってもFFTとして使えないか、…）、ここは素人考えで素数なら望ましい整数環になるんじゃなかろーか．．．などと勝手に決めつけました。
ではデカイ素数（値）をどうやってみつけるのか？ネットでミラー・ラビン素数判定法というのを見つけたので、十進BASICで試してみました。試しながらも、つらつらと考えました。この方法は被検査数値が素数でしたという確定結果を答えてくれません（たぶん素数でしょう、．．．みたいな。＝＝＞m(_ _)mすみません、決してこのアルゴリズムの批判とかでなく、素人の素朴な感想です。m(_ _)m）。それで他に簡単で確定結果を出してくれる方法は無いのか? 仕方なく、書店で探すことにしました。なるべく本の在庫が豊富そうな、池袋ジュンク堂書店へ。有りました !! 「[講座]数学の考え方⑯初等整数論」木田祐司（著）、この中に＜原始根テスト＞なる方法が書かれていました。
＜原始根テスト＞
或る数nが素数か？
＝＞n - 1を素因数分解する。 （例えば n - 1 = a1k1 * a2k2 * a3k3 だったとします）
＝＞その場合に或る数qを仮定して、nと各因数(a1、a2、a3)とqで作られた判定式をクリアしたら、
＝＞qはnの原始根（と言うらしい）で、nは素数である。
という感じです。ｑをどこまで探せば良いのか分かりませんが、時間がかかるので、2～1000に絞って、それで駄目なら検査終了（素数ではなさそう）とすることにしました。n －1の素因数分解は巨大な数では難しい訳ですが、別に暗号に使う為ではないので、n－1の因数が、判りやすく且つ種類が多くならないような数を候補として選びます。
これも十進BASICでササッと作成。探すのは、2(32＊i)より小さくて、出来るだけそれに近い素数。さらに都合の良い1の累乗根(回転子：複素数版のsin／cosに相当）が見つかるように、
2((32*i) - j) * k + 1、 2 <= k < 2j
の形の数値で上記原始根テストを行います。
素数が見つかったら回転子（整数環での1の累乗根）を求める訳ですが、数値がデカイですし何より累乗根を計算できないので、逆に数値を適当に選んでは累乗してみて1になるまでの限界を見極めて選びます。（全数の累乗はとても計算できないので、22、23、．．．、2xの各数値で累乗して、1になる限界をチェックします）。前述の様な素数が見つかった場合、大体2(32*i - j)乗で1になるω（回転子？）があるみたいです。
P=2u * v + 1
1≡ωw mod(P) (w=2u または w=2(u - 1))
というP、u、v、ωを探した訳です。その結果見つかった手頃なもので、整数環FFTを作って性能比較をしてみました。
64～320ビット同士の掛け算やモンゴメリ乗算は、大体c言語のインラインアセンブルコードを主体として作りました。
例１．P=228*13 + 1 (16進表示: 0xD0000001) (10進表示: 3489660929)
ω=11 (w=228)
1桁の単位=8ビット
例２．P=230*3 + 1 (16進表示: 0xC0000001) (10進表示: 3221225473)
ω=13 (w=230)
1桁の単位=8ビット
例３．P=259*27 + 1 (16進表示: 0xD80000000, 1)
(10進表示: 15564440312192434177)
ω=87 (w=259)
1桁の単位=16ビット
例４．P=292*7 + 1 (16進表示: 0x70000000, 0, 1)
(10進表示: 34662321099990647697175478273)
ω=11 (w=292)
1桁の単位=32ビット
例５．P=2121*81 + 1 (16進表示: 0xA2000000, 0, 0, 1)
(10進表示: 215334935317156371410416743765415821313)
ω=57 (w=2120)
1桁の単位=32ビット
例６．P=2189*3 + 1 (16進表示: 0x60000000, 0, 0, 0, 0, 1)
(10進表示: 2353913150770005286438421033702874906038383291674012942337)
ω=3 (w=2188)
1桁の単位=64ビット
例７．P=2316*13 + 1 (16進表示: 0xD0000000, 0, 0, 0, 0, 0, 0, 0, 0, 1)
(10進表示: 1735489466685739441945955136262761093114697424414780375581971306355553527196770446893656695635969)
ω=7 (w=2316)
1桁の単位=128ビット
性能比較：
Pが小さいと桁当たりの処理オーバヘッドが大きくなり、Pが大きいと手作り掛け算処理の速度が影響してきます。性能（と言ってもどんぐりの背比べです）的には、例４）P＝292*7+1、例６）P＝2189*3+1、例７）P＝2316 *13+1 がまあまあでした。浮動小数点（複素数）FFTといい勝負です。
KARATSUBAのところで書いたように、CPU（整数と浮動小数点の演算速度の違い）によってこの差も異なります。
今使ってるノートPCは、
（PC1）HP1000-1401（Win8、CPUはIntel Celeron CPU B830 CLK1.8GHz）
で、その前にメインで使っていたのは、
（PC2）ディスクトップＰＣ（ＰＣショップオリジナル）
（WinXP、CPUはAMD Athlon 64X2 Dual Core Proc 5400+ CLK2.8GHz）
ですが、PC2はPC1に比べて浮動小数点計算が相対的に遅いです。
いろいろ作った乗算処理（FFTとKARATSUBA方式）の速度（経過時間：秒）は、現在値としてこのくらいです。
我が家のPCでは、メモリ食いの浮動小数点（複素数）FFTを使って20million桁以上の計算は実行できません。
メモリ効率をみると、整数環FFTのほうが大勝?です。8ビットを処理するためにFFT処理で何ビット必要かを比べると、
標準（c言語の倍精度浮動小数点52ビット精度）： 128ビット（十進300万桁以上の計算の場合）
64ビット精度浮動小数点： 96ビット
手作り192ビット精度浮動小数点： 48ビット
整数環FFT（例４）： 24ビット
整数環FFT（例６）： 24ビット
整数環FFT（例７）： 20ビット
これから：
なにか改善?の気力が湧いて実行できたら、また記録を残したいと思います。
追記：
整数環FFT（一般のも含めて?）の回転子、入力桁の数を分離できるだけの値（桁数の累乗根）で良いのでしょうか?ヤル気が出たら調べてみましょう。しかし計算はちゃんとできてるようなので、今のところは問題ありません。この辺り、素人なので無駄なところに労力（神経）を割いてしまっているような気がします。
前の記事で、ラマヌジャン公式の計算ミスについて書きましたが、調べたら大事ではなく解決しました。最終式では”B”値を使っていない（チュドノフスキーの公式も同様）ので、最後の1項になるときの計算で、時間節約のためにB計算をしないようにしていたのですが、その判定条件を誤っていただけでした。よかった、よかった。
（注）本文は、数学素人の個人的な記録（単なるメモ）なので、誤りがある事などをご了解ください。
≪参考文献・記事≫
計算機代数入門（野呂正行）
数値計算の基礎（川上一郎）
「[講座]数学の考え方（16）初等整数論」木田祐司（著）朝倉書店
円周率.jp （ウィキペディア）
モンゴメリー乗算 （ウィキペディア）
FFT（高速フーリエ・コサイン・サイン変換）の概略と設計法
（Takuya OOURA）