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

「FFT を利用することで効率的に多倍長乗算の計算ができる」 ということはよく知られているがその詳細についての説明はあまり言及されない。 ここではその間をつなげる理論背景について説明する。

小難しい話にはなるが、多倍長整数の乗算に FFT を利用するためには まず多倍長整数をある程度の長さの数列に変換することになる。 Toom-Cook 法の導入部分で説明したような 多項式表現で考えてもらって構わない。 その上で DFT (FFT) を利用することで畳み込み乗算 (演算子として $*$ を使うことにする) が得られ、 再度多倍長整数の形へ変換することで計算が完了する。 \[ \begin{eqnarray} & A \mapsto \{a_j\},\quad B \mapsto \{b_j\}\\ & \{c_k\} = \{a_j\} * \{b_j\}\\ & \{c_k\} \mapsto C=A\times B \end{eqnarray} \]

前のページに書いたとおり FFT は DFT を効率的に計算するアルゴリズムで、 乗算との関係は DFT と畳み込み乗算の関係から導かれる。

Toom-Cook 法の導入部分で説明したような、 多項式同士の乗算をしたときのような数列の掛け合わせを 畳み込み乗算 (Convolution) という。 \[ \begin{eqnarray} C(x) &=& A(x)B(x)\\ &=& \left(\sum_{j=0}^{n-1}a_jx^j\right)\left(\sum_{j'=0}^{n-1}b_{j'}x^{j'}\right)\\ \sum_{k=0}^{2n-2} c_kx^k &=& \sum_{k=0}^{2n-2} \sum_{j+j'=k}a_jb_{j'}x^k\\ \end{eqnarray} \] \[ \therefore c_k = \sum_{j+j'=k}a_jb_{j'} \]

見ての通り、長さ $n$ の数列 $\{a_j\}$、$\{b_j\}$ の畳み込み乗算では長さ $2n-1$ の数列 $\{c_k\}$ が出来上がる。小学校でやった筆算のように書いてみると どのような計算がなされているか直感的に分かるだろう。 \[ \begin{array}{ccccc} &&&& a_{n-1} & a_{n-2} & \cdots & a_1 & a_0\\ &&& * & b_{n-1} & b_{n-2} & \cdots & b_1 & b_0\\ \hline &&&& a_{n-1}b_0 & a_{n-2}b_0 & \cdots & a_1b_0 & a_0b_0\\ &&& a_{n-1}b_1 & a_{n-2}b_1 & a_{n-3}b_1 & \cdots & a_0b_1 & \\ &&&& \vdots\\ & a_{n-1}b_{n-2} & \cdots & a_2b_{n-2} & a_1b_{n-2} & a_0b_{n-2}\\ a_{n-1}b_{n-1} & a_{n-2}b_{n-1} & \cdots & a_1b_{n-1} & a_0b_{n-1}\\ \hline c_{2n-2} & c_{2n-3} & \cdots & c_n & c_{n-1} & c_{n-2} & \cdots & c_1 & c_0 \end{array} \]

さてこのように $n$ 個のデータ同士の演算で $2n$ 個のデータ ($c_{2n-1}=0$ と考える) が生成されると計算理論的には群をなさないため色んな理論を適用することができない。 そこで多項式表現で言えば $\pmod{x^n-1}$ をすることで $n-1$ 次式 ($n$ 個のデータ) に落とし込む計算を考える。 \[ C(x) = A(x)B(x) \pmod{x^n-1} \] \[ c_k = \sum_{j=0}^{n-1}a_jb_{(k-j)\bmod n} = \sum_{j=0}^{k}a_jb_{k-j} + \sum_{j=k+1}^{n-1}a_jb_{k-j+n} \]

この計算を 巡回畳み込み乗算 (Cyclic Convolution) という。 (演算子には $\circledast$ や $\circledast_+$ を使うことにする) こちらも同様に筆算方式で書いてみると分かりやすい。 \[ \begin{array}{ccccc} & a_{n-1} & a_{n-2} & \cdots & a_1 & a_0\\ \circledast_+ & b_{n-1} & b_{n-2} & \cdots & b_1 & b_0\\ \hline & a_{n-1}b_0 & a_{n-2}b_0 & \cdots & a_1b_0 & a_0b_0\\ & a_{n-2}b_1 & a_{n-3}b_1 & \cdots & a_0b_1 & \color{red}{a_{n-1}b_1}\\ & \vdots & \vdots & & \color{red}\vdots & \color{red}\vdots\\ & a_0b_{n-1} & \color{red}{a_{n-1}b_{n-1}} & \color{red}\cdots & \color{red}{a_2b_{n-1}} & \color{red}{a_1b_{n-1}}\\ \hline & c_{n-1} & c_{n-2} & \cdots & c_1 & c_0 \end{array} \]

赤く書かれた部分は元の畳み込み乗算では $n\leqq k\lt 2n$ の $c_k$ に足しこまれていた部分をそのまま $n$ 個分右にずらして (多項式表現で言えば $n$ 次以上の係数を $n$ 次落として) 足し込んだ形になっている。

さらに一段階発展させて、巡回させた部分の符号を反転する 負巡回畳み込み乗算 (Nega-cyclic Convolution; 演算子は $\circledast_-$ を使うことにする) を紹介する。こちらは多項式表現で言えば $\pmod{x^n+1}$ での計算になる。 \[ C(x) = A(x)B(x) \pmod{x^n+1} \] \[ c_k = \sum_{j=0}^{k}a_jb_{k-j} - \sum_{j=k+1}^{n-1}a_jb_{k-j+n} \]

同様に筆算方式で書いてみると \[ \begin{array}{ccccc} & a_{n-1} & a_{n-2} & \cdots & a_1 & a_0\\ \circledast_- & b_{n-1} & b_{n-2} & \cdots & b_1 & b_0\\ \hline & a_{n-1}b_0 & a_{n-2}b_0 & \cdots & a_1b_0 & a_0b_0\\ & a_{n-2}b_1 & a_{n-3}b_1 & \cdots & a_0b_1 & \color{red}{-a_{n-1}b_1}\\ & \vdots & \vdots & & \color{red}\vdots & \color{red}\vdots\\ & a_0b_{n-1} & \color{red}{-a_{n-1}b_{n-1}} & \color{red}\cdots & \color{red}{-a_2b_{n-1}} & \color{red}{-a_1b_{n-1}}\\ \hline & c_{n-1} & c_{n-2} & \cdots & c_1 & c_0 \end{array} \]

となる。

改めて、DFT と畳み込み乗算の関係について説明しよう。 以下のように長さ $n$ の 2 数列 $\{a_j\}$、$\{b_j\}$ の DFT を考える。 \[ 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) \]

となる。 ここで ${\displaystyle \sum_{k=0}^{n-1}}\omega^{Jk}$ について \[ \sum_{k=0}^{n-1}\omega^{Jk}= \begin{cases} 0 & (J \bmod n \neq 0)\\ n & (J \bmod n = 0) \end{cases} \]

となるので \[ 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} \]

つまり巡回畳込み乗算を計算できたことになる。 しかし本来計算したかったものは巡回畳込み乗算ではなく 畳込み乗算である。 そこで元のデータのうち $n/2\leq j\lt n$ で $a_j=b_j=0$ であることを仮定して同じ計算(巡回畳込み乗算)を行うと以下のように 巡回する部分の値が 0 になり、 畳み込み乗算と同じ結果になることが分かる。

\[ \begin{array}{ccccc} & 0 & 0 & \cdots & 0 & a_{n/2-1} & a_{n/2-2} & \cdots & a_1 & a_0\\ \circledast_+ & 0 & 0 & \cdots & 0 & b_{n/2-1} & b_{n/2-2} & \cdots & b_1 & b_0\\ \hline & 0 & 0 & \cdots & 0 & a_{n/2-1}b_0 & a_{n/2-2}b_0 & \cdots & a_1b_0 & a_0b_0\\ & 0 & 0 & \cdots & a_{n/2-1}b_1 & a_{n/2-2}b_1 & a_{n/2-3}b_1 & \cdots & a_0b_1 & \color{red}{0}\\ & \vdots & \vdots & & \vdots & \vdots\\ & 0 & a_{n/2-1}b_{n/2-1} & \cdots & a_1b_{n/2-1} & a_0b_{n/2-1} & \color{red}{0} & \cdots & \color{red}{0} & \color{red}{0}\\ \hline & c_{n-1} & c_{n-2} & \cdots & c_{n/2} & c_{n/2-1} & c_{n/2-2} & \cdots & c_1 & c_0 \end{array} \]

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

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

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

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

DFT を使った計算の前に筆算形式で巡回畳み込み乗算の結果を導いておく。 \[ \begin{array}{rrrrrrrrr} (j: & 7 & 6 & 5 & 4 & 3 & 2 & 1 & 0)\\ & 31 & 41 & 59 & 26 & 53 & 58 & 97 & 93\\ \circledast_+ & 23 & 84 & 62 & 64 & 33 & 83 & 27 & 95\\ \hline & 2945 & 3895 & 5605 & 2470 & 5035 & 5510 & 9215 & 8835\\ & 1107 & 1593 & 702 & 1431 & 1566 & 2619 & 2511 & \color{red}{837}\\ & 4897 & 2158 & 4399 & 4814 & 8051 & 7719 & \color{red}{2573} & \color{red}{3403}\\ & 858 & 1749 & 1914 & 3201 & 3069 & \color{red}{1023} & \color{red}{1353} & \color{red}{1947}\\ & 3392 & 3712 & 6208 & 5952 & \color{red}{1984} & \color{red}{2624} & \color{red}{3776} & \color{red}{1664}\\ & 3596 & 6014 & 5766 & \color{red}{1922} & \color{red}{2542} & \color{red}{3658} & \color{red}{1612} & \color{red}{3286}\\ & 8148 & 7812 & \color{red}{2604} & \color{red}{3444} & \color{red}{4956} & \color{red}{2184} & \color{red}{4452} & \color{red}{4872}\\ & 2139 & \color{red}{713} & \color{red}{943} & \color{red}{1357} & \color{red}{598} & \color{red}{1219} & \color{red}{1334} & \color{red}{2231}\\ \hline & 27082 & 27646 & 28141 & 24591 & 27801 & 26556 & 26826 & 27075 \end{array} \]

次に 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} \mapsto A = \text{DFT}(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} \mapsto B = \text{DFT}(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 \mapsto = \text{DFT}^{-1}(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 \mapsto \{a_j\} = \{0, 0, 0, 0, 53, 58, 97, 93\}_{j=7..0} \] \[ 33832795 \mapsto \{b_j\} = \{0, 0, 0, 0, 33, 83, 27, 95\}_{j=7..0} \]

再度筆算で計算結果を求めておこう。 \[ \begin{array}{rrrrrrrrrl} (j: & 7 & 6 & 5 & 4 & 3 & 2 & 1 & 0)\\ & 0 & 0 & 0 & 0 & 53 & 58 & 97 & 93\\ \circledast_+ & 0 & 0 & 0 & 0 & 33 & 83 & 27 & 95\\ \hline & 0 & 0 & 0 & 0 & 5035 & 5510 & 9215 & 8835\\ & 0 & 0 & 0 & 1431 & 1566 & 2619 & 2511 & \color{red}{0}\\ & 0 & 0 & 4399 & 4814 & 8051 & 7719 & \color{red}{0} & \color{red}{0}\\ & 0 & 1749 & 1914 & 3201 & 3069 & \color{red}{0} & \color{red}{0} & \color{red}{0}\\ \hline & 0 & 1749 & 6313 & 9446 & 17721 & 15848 & 11726 & 8835 & \text{↓正規化}\\ & 18 & 13 & 09 & 24 & 80 & 66 & 14 & 35 \end{array} \]

今回は多倍長乗算との対応を確認するために畳み込み乗算の結果を出した後 1 要素 2 桁になるよう正規化しているが、その結果が求めたかった積 \[ 53589793 \times 33832795 = 1813092480661435 \]

と一致することが確認できる。

改めて 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} \mapsto A = \text{DFT}(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} \mapsto B = \text{DFT}(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 \mapsto c = \text{DFT}^{-1}(C) = \begin{pmatrix} 8835\\ 11726\\ 15848\\ 17721\\ 9446\\ 6313\\ 1749\\ 0 \end{pmatrix} \]

この結果は筆算の計算で見た、畳み込み乗算の結果と一致することが確認できる。