フーリエ変換
ここで唐突にフーリエ変換登場。フーリエ変換は（時間軸／空間軸に）リニアに並ぶ複素数列をフーリエ空間（周波数軸）に変換します。
二重for-loopで実現した総当たりの掛け算は、畳み込み（コンボリューション：comvolution）と呼ばれる数学的な操作です。2つのデータ列 a, b の畳み込みを a*b と表すことにします。
フーリエ変換には「畳み込み定理」と呼ばれる面白い性質があって、畳み込み a*b のフーリエ変換：F(a*b) は F(a) と F(b) の要素ごとの積:F(a)・F(b) と等しくなります。要素ごとの積ってことは a[0]*b[0], a[1]*b[1], ... なのですから、要素ごとの積に要する時間計算量はデータ数：N に比例、すなわちO(N)、畳み込みのO(N^2)に比べ圧倒的に小さくなります。
フーリエ空間上のデータ列を元の実空間に変換する逆フーリエ変換：F'(x)が存在しますから、畳み込み a*b をそのまま計算するんじゃなく:
- a と b のフーリエ変換：F(a), F(b) を求める
- F(a) と F(b) の要素ごとの積：F(a)・F(b) を求める
- F(a)・F(b) の逆フーリエ変換：F'(F(a)・F(b)) を求める
この3つのステップで F'(F(a)・F(b)) ⇒ F'(F(a*b)) ⇒ a*b が求まるってわけ。
フーリエ変換・逆フーリエ変換の時間計算量は高速フーリエ変換ならO(NlogN)です。
上記各ステップにおいて、フーリエ変換を2回、要素ごとの積を1回、逆フーリエ変換を1回行うので、つまるところ畳み込み演算 a*b を時間計算量:O(NlogN)で行える勘定になります。O(N^2)がO(NlogN)になるってことはバブルソートとクイックソートの速度比と同じ、Nが大きいほどその差は開いていきます。
……さてと、FFTによる巨大整数の掛け算を実装しましょうか。
コラム
今回使用するのはフーリエ変換ライブラリの定番「FFTW："Fastest Fourier Transform in the West（西側最速のフーリエ変換）"」、現時点での最新版はversion 3.3.4で、Windows用コンパイル済バイナリ・パッケージは32／64bitの両方が用意されています。
バイナリ・パッケージを展開するとヘッダとDLLがでてくるのですが、DLLのインポート・ライブラリ（LIB）が見当たりません。Visual C++で利用すべく、同梱されたDEFとDLLをlib.exeに食わせてLIBを作っておきましょう。コマンドラインから：
lib /machine:x86 /def:libfftw3-3.def lib /machine:x86 /def:libfftw3f-3.def
でライブラリ：libfftw3-3.lib, libfftw3f-3.libが生成されるので、ヘッダ（fftw3.h）とDLL, LIBを適当なディレクトリに配置しておいてください。2つのライブラリの違いは複素数の精度。実部/虚部の型がdoubleなのが前者、floatなのが後者です。64bitライブラリを生成する場合、libに与えるオプションを/machine:x64とします。
フーリエ変換による整数の掛け算：fftwf_multipliesのインターフェースは筆算方式のonpaper_multipliesと同じく、10進文字列std::string sa, sbを引数に与え、std::stringを返すことにしましょう。
std::string fftwf_multiplies(const std::string& sa, const std::string& sb) {
フーリエ変換の計算対象は複素数ですから、まず与えられた10進文字列を各桁ごとの複素数列：vector<complex<float>>に変換します。N桁の2数の積は最大2N桁となるので、複素数列の長さは与えられた2数の桁数の大きいほうの2倍とします。ついでに結果を納める複素数列も同じサイズで用意します。
using real_type = float;
using complex_type = complex<real_type>;
int digits = (int)max(sa.size(), sb.size()) * 2;
// ca, cb, cc : a, b, c の 複素数表現
complex_type* ca = (complex_type*)fftwf_alloc_complex(digits);
complex_type* cb = (complex_type*)fftwf_alloc_complex(digits);
complex_type* cc = (complex_type*)fftwf_alloc_complex(digits);
// 文字列sa,sbを（下位の桁から順に）複素数列ca,cbに変換（虚部は0,余った桁は 0+0i で埋める）
auto char2complex = [](char ch) { return complex_type((real_type)(ch - '0'));};
fill(transform(rbegin(sa), rend(sa), ca, char2complex), ca+digits, complex_type());
fill(transform(rbegin(sb), rend(sb), cb, char2complex), cb+digits, complex_type());
前述のステップに従い、caとcbにフーリエ変換を施し、その結果の要素ごとの積をccに求め、さらにccに対し逆フーリエ変換をかけます。
fftwf_plan plan_a = fftwf_plan_dft_1d(digits, (fftwf_complex*)ca, (fftwf_complex*)ca,
FFTW_FORWARD, FFTW_ESTIMATE); // フーリエ変換
fftwf_plan plan_b = fftwf_plan_dft_1d(digits, (fftwf_complex*)cb, (fftwf_complex*)cb,
FFTW_FORWARD, FFTW_ESTIMATE); // フーリエ変換
fftwf_plan plan_c = fftwf_plan_dft_1d(digits, (fftwf_complex*)cc, (fftwf_complex*)cc,
FFTW_BACKWARD, FFTW_ESTIMATE); // 逆フーリエ変換
// ca, cb をフーリエ変換
fftwf_execute(plan_a);
fftwf_execute(plan_b);
// cc = ca * cb
transform(ca, ca+digits, cb, cc, multiplies<complex_type>());
// cc を逆フーリエ変換
fftwf_execute(plan_c);
FFTWによる逆フーリエ変換の結果は正規化されていないので、各要素の値を要素数で割って正規化し、さらに複素数列から整数列に変換します。
// 実数部を桁数(digits)で割って正規化
vector<int> c(digits,0);
transform(cc, cc+digits, begin(c),
[=](complex_type z) { return (int)(z.real()/digits+0.5f);});
あとは筆算と同様、繰り上げ処理を行い、文字列化して完了。
// 桁あふれの繰り上げ
int carry = 0;
transform(begin(c), end(c), begin(c),
[&](int n) { int t = n + carry; carry = t /10; return t % 10; });
// 文字列化(ゼロ・サプレスを行う)
string result;
result.reserve(digits);
bool nonzero = false;
for_each(rbegin(c), rend(c),
[&](int ch) { if ( ch != 0 && !nonzero ) { nonzero = true; }
if ( nonzero ) { result.push_back(ch + '0'); } });
return result;
筆算版とフーリエ変換版それぞれの畳み込み演算部に所要時間の計測ルーチンを埋め込み、いくつかの桁数で掛け算を行って所要時間を計ってみました。csvで出力しexcelに食わせて描いたグラフがコレ（フーリエ変換版は単精度（complex<float>）と倍精度（complex<double>）の2つを用意しました）。
単精度だと演算精度が足りず、桁数の大きな掛け算で筆算版とは異なった結果になっちゃいました（所要時間：0はそのシルシ）。60桁を超えたあたりで筆算とフーリエ変換との所要時間が逆転し、その後その差は開く一方です。ソートと同様、アルゴリズムによるパフォーマンスの差を体感できる良い例となりました。
