DFTと畳み込み乗算と多倍長乗算

長さ $n$ の 2 数列 $\{x_i\}$ と $\{y_i\}$ があるとき、その畳み込み乗算 $\{z_i\}$

z_k = \sum_{i+j \equiv k \pmod{n}} x_i y_j

は離散フーリエ変換(DFT) を使って

\{z_k\} = {\rm DFT}^{-1} ({\rm DFT}(\{x_i\}) * {\rm DFT}(\{y_j\}))

と計算できる。 $*$ については同じインデックスの要素同士をかける "項別乗算" である。概念的に C++ コードで表すと以下のようになる。

// x と y の畳み込み演算結果を z に格納する
void Convolution(Element x[], Element y[], int64 n, Element z[]) {
  DFT(x, n, Forward, xf);  // xf = DFT(x)
  DFT(y, n, Forward, yf);  // yf = DFT(x)
  for (int64 i = 0; i < n; ++i)
    zf[i] = xf[i] * yf[i];
  DFT(zf, n, Backward, z); // z = DFT^-1(zf)
}

コンピュータでの計算については、特に DFT の演算に FFT アルゴリズムを使うことで計算量が $O(n\log n)$ になること、畳み込み演算結果を正規化すれば多倍長整数の乗算となること、から多倍長整数の乗算でよく用いられている。

Integer& operator*(const Integer& x, const Integer& y) {
  const int64 n = ...;  // 適切に決める。常に定数でも可。
  Element xx[n], yy[n];
  SplitIntoElements(x, n, xx);  // 区切って数列に分ける。データ構造次第では必要ない
  SplitIntoElements(y, n, yy);
  Element zz[n];
  Convolution(xx, yy, n, zz);
  Integer z(zz, n);  // 正規化(後述) & 多倍長整数へのデータ変換。
  return z;
}

FFT の基本原理

まずはざっくりと FFT が $O(n \log n)$ の計算量で行える原理を記しておく。数式としては

X_k = \sum_{j=0}^{n-1} x_j w^{jk}
\quad \left(x_j = n^{-1} \sum_{k=0}^{n-1} X_k w^{-jk}\right)

の定義を

\begin{eqnarray}
X_k &= \sum_{j=0}^{n/2-1} x_{2j} w^{2jk} + w^k\sum_{j=0}^{n/2-1} x_{2j+1} w^{2jk}\\
X_{k+n/2} &= \sum_{j=0}^{n/2-1} x_{2j} w^{2jk} + w^{n/2}w^k\sum_{j=0}^{n/2-1} x_{2j+1} w^{2jk} 
\end{eqnarray}

と変形することで、サイズ $n/2$ の DFT ２回分に分割する。これによりサイズ $n$ の計算量 $F(n)$ は $F(n)=2F(n/2)+\alpha n$ と漸化式で表わすことができ、再帰展開することで $F(n)=O(n\log n)$ となる。ちなみに、後半の式では $w^{n/2}=-1$ ということを利用することで更にシンプルに

\begin{eqnarray}
X_k &= \sum_{j=0}^{n/2-1} x_{2j} w^{2jk} + w^k\sum_{j=0}^{n/2-1} x_{2j+1} w^{2jk}\\
X_{k+n/2} &= \sum_{j=0}^{n/2-1} x_{2j} w^{2jk} - w^k\sum_{j=0}^{n/2-1} x_{2j+1} w^{2jk} 
\end{eqnarray}

とすることができ、一般にはこの式ばかりが世の中で知れ渡っているが計算量を落とす原理としては前の分割式で問題ない。(というよりも2基底以外の FFT を構成するためには前の式の形が重要なのだが、別の話なので今回は置いておく。)

さて、このような演算量の落とし方をこの記事の中では "FFTアルゴリズム" ということにしよう。すると、FFT アルゴリズムを適用するためには以下に記したルールが守られていればいいことがわかる。

原始 $n$ 乗根 $w$ が存在すること。これは数学的には $w^n=1$、 $w^k \neq 1 (0 < k < n)$ ということである。
必須ではないが、$w^{(n/2)}=-1$ であると上記の理由で計算そのものを簡略化することができる。
$w^{-1}$、$n^{-1}$ が簡単に計算できること。

多倍長乗算における基本 : DFT を利用した畳み込み乗算の制限

従来の基本原理では $O(n\log n)$ の効率になり、筆算方式の $O(n^2)$ からすれば効率的になっているが、コンピュータでの実装を見ると以下のようなデメリットがある。

入力データが実数なので複素数の虚数部分で無駄にメモリを食う
浮動小数点数を使うので、厳密な計算ができず、整数丸めのためにさらにビットを必要とする
倍精度浮動小数点数 (double) の仮数部は実質 53 bit の有効ビットで、11 bit の指数部が無駄になっている

2, 3 点目について簡単な実験で整数化する際の誤差を計測してみたところ、以下の表のようになった。実験の詳細としては

1 要素 (double) に 16 bit ずつのデータをランダムに配置。
畳み込み乗算結果が $n$ 要素分になるよう、元のデータの [n/2 .. n-1] は 0
各 $n$ について 10 回ずつ演算をし、"整数化誤差の最大値" の平均を見る

$n$	平均誤差
4096	0.000e+00
8192	8.789e-04
16384	1.563e-03
32768	7.617e-03
65536	1.094e-02
131072	3.438e-02
262144	7.031e-02
524288	1.625e-01
1048576	3.688e-01
2097152	5.000e-01

4n 乗根を用いた DWT

ここで紹介する方法は 1 つめのデメリットに対する解決策である。簡単に言えば複素数の虚部にも使えるデータを入れることで無駄をなくす手法である。

$\{x_i\}$ と $\{y_i\}$ の長さの合計が $2n$ 以下となるように $n$ を設定する。
両者とも後ろに 0 を埋めて $2n$ の長さにする。
以下のようにまとめて長さ $n$ の複素数列 $\{x'_j\}$、$\{y'_j\}$ を定義する。

x'_j = (x_j + i x_{j+n})w^{j/4},\quad y'_j = (y_j + i y_{j+n})w^{j/4}

ちなみに $w$ は 1 の原始 $n$ 乗根である。するとこれらの畳み込み乗算

\{z'_k\} = {\rm DFT}^{-1} ({\rm DFT}(\{x'_i\}) * {\rm DFT}(\{y'_j\}))

は同様に

z'_j = (z_j + i z_{j+n})w^{j/4},\quad z_j + i z_{j+n}=z'_j w^{-j/4}

という形になっているので、必要なメモリ量は半分に、演算量もほぼ半分＋α ($w^{j/4}$ をかける処理、データの並べ替え)になる。

実部と虚部をまとめた DFT

先ほどの例と同様に複素数の虚部にもデータを置くことで容量の無駄を省く方法であるが、データの並べ方を以下のようにする。

x'_j = x_{2j} + i x_{2j+1},\quad y'_j = y_{2j} + i y_{2j+1}

このとき、 $\{X_k\} = {\rm DFT}(\{x_i\})$ と $\{X'_k\} = {\rm DFT}(\{x'_i\})$ との関係が

\begin{cases}
A'_k = \dfrac{1}{2}(X'_k - \overline{X'_{n/2-k}})(1 + iw^k)\\
X_k = X'_k - A'_k, \quad X_{n/2-k} = X'_{n/2-k} + \overline{A'_k}\\
X_0 = \Re(X'_0) + \Im(X'_0), \quad X_{n/2}=\Re(X'_0)-\Im(X'_0)
\end{cases}

となる。ちなみに逆方向の変換は

\begin{cases}
A_k = \dfrac{1}{2}(X_k-\overline{X_{n/2-k}})(1-iw^{-k})\\
X'_k = X_k - A_k, \quad X'_{n/2-k}=X_{n/2-k}+\overline{A_k}\\
X'_0=\dfrac{1}{2}(X_0 + X_{n/2}) + \dfrac{i}{2}(X_0 - X_{n/2})
\end{cases}

である。この方法では一般的なデータの並び順(Little endian)をそのまま複素数として取り扱うことができるので、C/C++ でいう union のような形で捉えればデータの並び替え処理が必要ない。

整数環で行う DFT

整数の範囲でも FFT アルゴリズムを適用できる計算があり、 1 つめ、 2 つめのデメリットは完全になくなる。単純に言えば元の条件を満たす体や環で計算を行えばいいだけである。

p < 2^64 で例

コンピュータでの計算メリットがある例として $\bmod, p=2^{64}-2^{32}+1$ の計算を挙げる。 (参考;英語)

この環では 7 が原始 $p-1$ 乗根であり、$p-1=2^{32} \cdot 3 \cdot 5 \cdot 17 \cdot 257 \cdot 65537$ ということから、 $n=2^{32}$ までの 2 基底 FFT に加え、３基底や 5 基底を含めることができる。

加えて 2 はこの環で原始 192 乗根 (ただし $7^{(p-1)/192} \bmod p\not=2$) となっているので、

\begin{cases}
2^{64} \equiv 2^{32} - 1 \pmod p\\
2^{96} \equiv -1 \pmod p\\
2^{128} \equiv -2^{32} \pmod p\\
2^{192} \equiv 1 \pmod p\\
2^n \cdot 2^{192-n} \equiv 1 \pmod p
\end{cases}

となることから、 64bit の数を 32bit ずつに分割して計算すると $\pmod p$ の計算が楽にできる。

\begin{eqnarray}
&& 2^{96}x_3 + 2^{64}x_2 + 2^{32}x_1 + x_0 \bmod p\\
&=& 2^{32}(x_1+x_2) + x_0 - x_3 - x_2
\end{eqnarray}

Schönhage-Strassen

3 つめのデメリットは、要するに畳み込み乗算結果が最大で

(要素数)(基数 - 1)^2

となるため、それをレジスタ単位で格納するのに限界があるということである。つまり解決する一つの素直な方法は各要素の計算にも多倍長計算を導入することであり、Schönhage-Strassen の方法はそれを解決する一つの方法である。

細かい計算方法については別記事に纏めているので、そちらを参照してほしい。