円周率.jp (http://xn--w6q13e505b.jp/method/fft/2dfft.html)

2次元化FFT

多数桁乗算に用いられる FFT で計算対象になる数値は, 一般に桁数 $n$ の 1 次元配列でもって表される. FFT の演算ではほぼランダムアクセスに近いメモリアクセスをするため, $n$ がキャッシュメモリより大きくなると極端に演算速度が遅くなる.

これを回避するため,$n=qp$ と因数分解できると仮定し,FFT の計算式を次のように書き換える.

\[ k = k_1p + k_2,\ j = j_1q + j_2\quad (0\leq k_1, j_2 \lt q,\ 0\leq k_2, j_1 \lt p) \] \[ \begin{eqnarray} A_{k_1p+k_2} &=& \sum_{j_2=0}^{q-1}\sum_{j_1=0}^{p-1} a_{j_1q+j_2}\omega^{(j_1q+j_2)(k_1p+k_2)}\\ &=& \sum_{j_2=0}^{q-1}\sum_{j_1=0}^{p-1} a_{j_1q+j_2}\omega^{j_1k_2q+j_2k_1p+j_2k_2}\\ &=& \sum_{j_2=0}^{q-1} \left( \left(\sum_{j_1=0}^{p-1} a_{j_1q+j_2}\omega_p^{j_1k_2} \right) \omega^{j_2k_2} \right) \omega_q^{j_2k_1}\\ \end{eqnarray} \]

このように変形することで,$p$,$q$ 程度の長さの FFT を繰り返すことになるので,$p$,$q$ 次第ではキャッシュメモリに収まり演算全体の高速化が期待できる. 基本的には max ($p,q$) を最小化するため, $p\simeq q \simeq \sqrt{n}$ にする.

2 次元化 (four-step FFT / six-step FFT)

上記の式では若干わかりにくいので,入力となる 1 次元データ列を 2 次元に配置し再解釈することで分かりやすくなる. 大きく 4 つの計算ステップを踏むことから four-step FFT という.

$a_0[j_1,j_2] = a_{j_1q+j_2}$ (配列の 2 次元再解釈)
1st STEP $a_1[k_2,j_2] = {\displaystyle \sum_{j_1=0}^{p-1}} a_0[j_1,j_2]\omega_p^{j_1k_2}$ $({\rm for}\ 0\leq j_2\lt q)$
2nd STEP $a_2[k_2, j_2] = a_1[k_2, j_2] \omega^{j_2k_2}$ $({\rm for}\ 0\leq k_2\lt p, 0\leq j_2 \lt q)$
3rd STEP $a_3[j_2, k_2] = a_2[k_2, j_2]$ $({\rm for}\ 0\leq k_2\lt p, 0\leq j_2 \lt q)$
4th STEP $a_4[k_1,k_2] = {\displaystyle \sum_{j_2=0}^{q-1}}a_3[j_2, k_2] \omega_q^{j_2k_1}$ $({\rm for}\ 0\leq k_2\lt p)$
$a_{k_1p+k_2} = a_4[k_1,k_2]$ (配列の 1 次元再解釈)

four-step FFT を素直に解釈すると各 FFT を行う際に用いるデータが不連続に並んでいるため, 計算効率が非常に悪く見える. それを改善するために,最初と最後(厳密には最後でなくても良い)に行列を転置することで FFT でアクセスするメモリが連続になるようにしたものを six-step FFT という.

3 次元化 FFT (nine-step FFT)

同様に 3 次元化も可能である.2 次元化で使う 1 次元 FFT がキャッシュサイズに載らない場合にはこちらを使うことで効果が上がるかも知れない. また,分散メモリ型並列計算にも応用させやすい. 式変形は面倒だが,1 度やっておくと理解が深まる. 変換前後のデータ $a_j$,$A_k$ の 3 軸をどう取るかについては, 変換後の式を考えて $\omega$ の項が少なくなるように取ると良い. 1 つの分割,計算例を以下に書く.

\[ n = pqr, k = (k_1q + k_2)p + k_3, j = (j_1q + j_2)r + j_3 \] \[ (0\leq k_3,j_1\lt p,\ 0\leq k_2,j_2 \lt q,\ 0\leq k_1,j_3\lt r) \] \[ \begin{eqnarray} A_{k_1pq+k_2p+k_3} &=& \sum_{j_1=0}^{p-1}\sum_{j_2=0}^{q-1}\sum_{j_3=0}^{r-1} a_{j_1qr+j_2r+j_3} \omega^{(j_1qr+j_2r+j_3)(k_1pq+k_2p+k_3)}\\ &=& \sum_{j_1=0}^{p-1}\sum_{j_2=0}^{q-1}\sum_{j_3=0}^{r-1} a_{j_1qr+j_2r+j_3} \omega^{j_1k_3qr+j_2k_2pr+j_3k_1pq+j_2k_3r+j_3k_2p+j_3k_3}\\ &=& \sum_{j_3=0}^{r-1}\left(\left(\sum_{j_2=0}^{q-1}\left(\left(\sum_{j_1=0}^{p-1} a_{j_1qr+j_2r+j_3} \omega_p^{j_1k_3} \right) \omega^{(j_2r+j_3)k_3}\right) \omega_q^{j_2k_2}\right)\omega_{qr}^{j_3k_2}\right)\omega_r^{j_3k_1} \end{eqnarray} \]

行列転置を含めて再解釈

3 次元化 FFT では最初と最後を含めて 4 回の行列転置が必要になる. また,分散メモリ型計算機を用いてこの計算を行う場合はそのうち 何回かはデータ通信を用いての転置となるため,通信速度がネックになりやすい.

$a_0[j_1,j_2,j_3]=a_{j_1qr+j_2r+j_3}$ (配列の 3 次元再解釈)
(1) $a_1[j_3,j_2,j_1] = a_0[j_1,j_2,j_3]$ $({\rm for}\ 0\leq j_1\lt p,\ 0\leq j_2\lt q,\ 0\leq j_3\lt r)$
(2) $a_2[j_3,j_2,k_3]={\displaystyle \sum_{j_1=0}^{p-1}} a_1[j_3,j_2,j_1]\omega_p^{j_1k_3}$ $({\rm for}\ 0\leq j_2\lt q,\ 0\leq j_3\lt r)$
(3) $a_3[j_3,j_2,k_3] = a_2[j_3,j_2,k_3]\omega^{(j_2r+j_3)k_3}$ $({\rm for}\ 0\leq k_3\lt p,\ 0\leq j_2\lt q,\ 0\leq j_3\lt r)$
(4) $a_4[j_3,k_3,j_2] = a_3[j_3,j_2,k_3]$ $({\rm for}\ 0\leq k_3\lt p,\ 0\leq j_2\lt q,\ 0\leq j_3\lt r)$
(5) $a_5[j_3,k_3,k_2] = {\displaystyle \sum_{j_2=0}^{q-1}} a_4[j_3,k_3,j_2] \omega_q^{j_2k_2}$ $({\rm for}\ 0\leq k_3\lt p,\ 0\leq j_3\lt r)$
(6) $a_6[j_3,k_3,k_2] = a_5[j_3,k_3,k_2] \omega_{qr}^{j_3k_2}$ $({\rm for}\ 0\leq k_3\lt p,\ 0\leq k_2\lt q,\ 0\leq j_3\lt r)$
(7) $a_7[k_3,k_2,j_3] = a_6[j_3,k_3,k_2]$ $({\rm for}\ 0\leq k_3\lt p,\ 0\leq k_2\lt q,\ 0\leq j_3\lt r)$
(8) $a_8[k_3,k_2,k_1] = {\displaystyle \sum_{j_3=0}^{r-1}} a_7[k_3,k_2,j_3] \omega_r^{j_3k_1}$ $({\rm for}\ 0\leq k_3\lt p,\ 0\leq k_2\lt q)$
(9) $a_9[k_1,k_2,k_3] = a_8[k_3,k_2,k_1]$ $({\rm for}\ 0\leq k_3\lt p,\ 0\leq k_2\lt q,\ 0\leq k_1\lt r)$
$a_{k_1pq+k_2p+k_3} = a_9[k_1,k_2,k_3]$ (配列の 1 次元再解釈)