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

フーリエ変換と乗算

FFT(DFT) を利用した乗算についての基礎を述べる. 2 つの多倍長数 $\{a_j\}$,$\{b_j\}$ を DFT した結果をそれぞれ $\{A_k\}$,$\{B_k\}$ とする.

\[ A_k = \sum_{j=0}^{n-1} a_j \omega^{jk},\quad B_k = \sum_{j=0}^{n-1} b_j \omega^{jk} \]

これらを項別に掛け合わせたデータ列を $\{C_k\}$ とする.

\[ C_k = A_kB_k = \left(\sum_{j=0}^{n-1}a_j\omega^{jk}\right) \left(\sum_{j'=0}^{n-1}b_j\omega^{j'k}\right) \]

この$\{C_k\}$ に IDFT をかけて出てくるデータ列 $\{c_j\}$ を計算すると

\[ c_l = \frac{1}{n}\sum_{k=0}^{n-1} C_k\omega^{-lk} = \frac{1}{n}\sum_{k=0}^{n-1}\left( \left(\sum_{j=0}^{n-1}a_j\omega^{jk}\right) \left(\sum_{j'=0}^{n-1}b_j\omega^{j'k}\right) \omega^{-lk}\right) \]

ここで,$J\not\equiv0 \pmod{n}$ のとき ${\displaystyle \sum_{k=0}^{n-1}}\omega^{Jk}=0$ で,$J\equiv0 \pmod{n}$ のときは ${\displaystyle \sum_{k=0}^{n-1}}\omega^{Jk}=n$ となるので

\[ c_l = \sum_{j=0}^{n-1} a_jb_{(l-j)\bmod n} = \sum_{j=0}^{l} a_jb_{l-j} + \sum_{j=l+1}^{n-1} a_jb_{l-j+n} \tag{1} \]

この式はインデックスを求めてみるとわかる通り図 1 の様な演算をしており,これを「(巡回)畳込み乗算」という. 図中の緑色の部分は ${\displaystyle \sum_{j=0}^{l}}$ の部分で, 黄色の部分は ${\displaystyle \sum_{j=l+1}^{n-1}}$ で計算される.

図 9: 畳込み乗算

このうち,$n/2\leq j\lt n$ で $a_j=b_j=0$ であれば図 2 のようになるので, これはつまり掛け算の結果と同様になることを示している.

図 10: 上半分が0の畳込み乗算

具体例

文字式だけでは分かりにくいので具体的な計算を,途中経過を含めて記述する.

全要素に値を入れている場合

まずは全要素に値を入れ,結果が畳込み乗算になることを確認する. 計算対象となる式を決め,計算の準備として値を 1 要素に 2 桁ずつ詰め込み, 多倍長数の形式に分割する.

\[ 3141592653589793 \times 2384626433832795 \] \[ 3141592653589793 \rightarrow \{a_j\} = \{31, 41, 59, 26, 53, 58, 97, 93\}_{j=7..0} \] \[ 2384626433832795 \rightarrow \{b_j\} = \{23, 84, 62, 64, 33, 83, 27, 95\}_{j=7..0} \]

DFT を使った計算の前に,筆算形式で畳込み乗算の結果を導いておく. 一般的な乗算ではこの枠より左側に伸びる部分を, 回り込ませて黄色の部分へ移動させる. その結果,最後に求める和が単純な乗算結果と異なることとなる. ちなみに,$3141592653589793 \times 2384626433832795 = 7491524886085135380744680661435$ である.

k 7654 3210
31415926 53589793
(×) 23846264 33832795
 
29453895 56052470 50355510 92158835
11071593 7021431 15662619 2511 837
48972158 43994814 80517719 25733403
8581749 19143201 30691023 13531947
33923712 62085952 19842624 37761664
35966014 57661922 25423658 16123286
81487812 26043444 49562184 44524872
2139 713 9431357 5981219 13342231
 
27082276462814124591 27801265562682627075

次に DFT を使った計算の過程を追っていく. $n=8$ なので $\omega=\dfrac{1}{\sqrt{2}}-\dfrac{i}{\sqrt{2}}$ である。まずは $\{a_j\}$,$\{b_j\}$ の DFT を求める.

\[ a = \begin{pmatrix} 93\\ 97\\ 58\\ 53\\ 26\\ 59\\ 41\\ 31 \end{pmatrix} \rightarrow A = \begin{pmatrix} 458\\ 67+\frac{16}{\sqrt2}-(17+\frac{60}{\sqrt2})i\\ 20-72i\\ 67-\frac{16}{\sqrt2}+(17-\frac{60}{\sqrt2})i\\ -22\\ 67-\frac{16}{\sqrt2}-(17-\frac{60}{\sqrt2})i\\ 20+72i\\ 67+\frac{16}{\sqrt2}+(17+\frac{60}{\sqrt2})i \end{pmatrix} \] \[ b = \begin{pmatrix} 95\\ 27\\ 83\\ 33\\ 64\\ 62\\ 84\\ 23 \end{pmatrix} \rightarrow B = \begin{pmatrix} 471\\ 31-\frac{45}{\sqrt2}+(1+\frac{25}{\sqrt2})i\\ -8-33i\\ 31+\frac{45}{\sqrt2}-(1-\frac{25}{\sqrt2})i\\ 181\\ 31+\frac{45}{\sqrt2}+(1-\frac{25}{\sqrt2})i\\ -8+33i\\ 31-\frac{45}{\sqrt2}-(1+\frac{25}{\sqrt2})i \end{pmatrix} \]

次に要素ごとに積を取り $\{C_k\}$ を計算する.

\[ C = \begin{pmatrix} 458 \cdot 471\\ \left(67+\frac{16}{\sqrt2}-(17+\frac{60}{\sqrt2})i\right) \left(31-\frac{45}{\sqrt2}+(1+\frac{25}{\sqrt2})i\right)\\ (20-72i)(-8-33i)\\ \left(67-\frac{16}{\sqrt2}+(17-\frac{60}{\sqrt2})i\right) \left(31+\frac{45}{\sqrt2}-(1-\frac{25}{\sqrt2})i\right)\\ -22 \cdot 181\\ \left(67-\frac{16}{\sqrt2}-(17-\frac{60}{\sqrt2})i\right) \left(31+\frac{45}{\sqrt2}+(1-\frac{25}{\sqrt2})i\right)\\ (20+72i)(-8+33i)\\ \left(67+\frac{16}{\sqrt2}+(17+\frac{60}{\sqrt2})i\right) \left(31-\frac{45}{\sqrt2}-(1+\frac{25}{\sqrt2})i\right) \end{pmatrix} = \begin{pmatrix} 215718\\ 2484-\frac{2034}{\sqrt2} + (1090+\frac{596}{\sqrt2})i\\ -2536 - 84i\\ 2484+\frac{2034}{\sqrt2} - (1090-\frac{596}{\sqrt2})i\\ -3982\\ 2484+\frac{2034}{\sqrt2} + (1090-\frac{596}{\sqrt2})i\\ -2536 + 84i\\ 2484-\frac{2034}{\sqrt2} - (1090+\frac{596}{\sqrt2})i \end{pmatrix} \]

最後に IDFT をかけて $\{c_j\}$ を求める.

\[ C \rightarrow c = \begin{pmatrix} 27075\\ 26826\\ 26556\\ 27801\\ 24591\\ 28141\\ 27646\\ 27082 \end{pmatrix} \]

この結果は確かに畳込み乗算の結果と一致することが確認できる.

上半分が 0 の場合

次に,基本的な流れを同じくしたまま,$a_{4..7}$,$b_{4..7}$ に 0 を入れて乗算結果と比較する. 計算する前に,計算式と計算に使う配列要素を確認しておく.

\[ 53589793 \times 33832795 \] \[ 53589793 \rightarrow \{a_j\} = \{0, 0, 0, 0, 53, 58, 97, 93\}_{j=7..0} \] \[ 33832795 \rightarrow \{b_j\} = \{0, 0, 0, 0, 33, 83, 27, 95\}_{j=7..0} \]

また,筆算で計算結果を求めておく.ここでは上記畳込み乗算との対比を明確にするため 0 が代入されている $b_{4..7}$ との積も明記しておく. 最後に和を取る部分では先に要素ごとの和を求め,その後で 1 要素 2 桁になるよう正規化している.

$j$ : 7654 3210
0000 53589793
$\times$ 0000 33832795
 
00 00 50355510 92158835
00 01431 15662619 25110
00 43994814 80517719 00
01749 19143201 30690 00
00 00 00 00
00 00 00 00
00 00 00 00
00 00 00 00
 
0174963139446 1772115848117268835 ↓正規化
18130924 80661435

次に DFT を使った計算の過程を追っていく. こちらも $n=8$ なので $\omega=\dfrac{1}{\sqrt2}-\dfrac{1}{\sqrt2}i$ となる. まずは $a_j$,$b_j$ の DFT を求める.

\[ a = \begin{pmatrix} 93\\ 97\\ 58\\ 53\\ 0\\ 0\\ 0\\ 0 \end{pmatrix} \rightarrow A = \begin{pmatrix} 338\\ 93+\frac{44}{\sqrt2}-(95+\frac{150}{\sqrt2})i\\ -2-44i\\ 93-\frac{44}{\sqrt2}+(95-\frac{60}{\sqrt2})i\\ 38\\ 93-\frac{44}{\sqrt2}-(95-\frac{60}{\sqrt2})i\\ -2+44i\\ 93+\frac{44}{\sqrt2}+(95+\frac{150}{\sqrt2})i\\ \end{pmatrix} \] \[ b = \begin{pmatrix} 95\\ 27\\ 83\\ 33\\ 0\\ 0\\ 0\\ 0 \end{pmatrix} \rightarrow B = \begin{pmatrix} 238\\ 95-\frac{6}{\sqrt2}-(83+\frac{60}{\sqrt2})i\\ 12+6i\\ 95+\frac{6}{\sqrt2}+(83-\frac{60}{\sqrt2})i\\ 118\\ 95+\frac{6}{\sqrt2}-(83-\frac{60}{\sqrt2})i\\ 12-6i\\ 95-\frac{6}{\sqrt2}+(83+\frac{60}{\sqrt2})i\\ \end{pmatrix} \]

次に要素ごとに積を取り $\{C_k\}$ を計算する.

\[ C = \begin{pmatrix} 338 \cdot 238\\ \left(93+\frac{44}{\sqrt2}-(95+\frac{150}{\sqrt2})i\right) \left(95-\frac{6}{\sqrt2}-(83+\frac{60}{\sqrt2})i\right)\\ (-2-44i)(12+6i)\\ \left(93-\frac{44}{\sqrt2}+(95-\frac{60}{\sqrt2})i\right) \left(95+\frac{6}{\sqrt2}+(83-\frac{60}{\sqrt2})i\right)\\ 38 \cdot 118\\ \left(93-\frac{44}{\sqrt2}-(95-\frac{60}{\sqrt2})i\right) \left(95+\frac{6}{\sqrt2}-(83-\frac{60}{\sqrt2})i\right)\\ (-2+44i)(12-6i)\\ \left(93+\frac{44}{\sqrt2}+(95+\frac{150}{\sqrt2})i\right) \left(95-\frac{6}{\sqrt2}+(83+\frac{60}{\sqrt2})i\right)\\ \end{pmatrix} = \begin{pmatrix} 80444\\ -3682-\frac{14528}{\sqrt2} - (17614+\frac{22912}{\sqrt2})i\\ 240-540i\\ -3682+\frac{14528}{\sqrt2} + (17614-\frac{22912}{\sqrt2})i\\ 4484\\ -3682+\frac{14528}{\sqrt2} - (17614-\frac{22912}{\sqrt2})i\\ 240+540i\\ -3682-\frac{14528}{\sqrt2} + (17614+\frac{22912}{\sqrt2})i\\ \end{pmatrix} \]

最後に,IDFT をかけて $\{c_j\}$ を求める.

\[ C \rightarrow c = \begin{pmatrix} 8835\\ 11726\\ 15848\\ 17721\\ 9446\\ 6313\\ 1749\\ 0 \end{pmatrix} \]

この結果は筆算の,要素ごとの和を出した段階と一致することが確認できる.