math
数学
FFT
多倍長整数

FFT を用いた多倍長整数の乗算の効率化 (2)

More than 1 year has passed since last update.

この記事はFFT を用いた多倍長整数の乗算の効率化 (1) からの派生記事なので、それを読んでいる前提に書いてます。


畳み込み乗算

多倍長乗算と DFT との架け橋として使われる畳み込み乗算は、詳しくは巡回畳み込み乗算と言われるもので、準備している要素数を超えた部分の乗算結果が下位に巡回して足し込まれる 正巡回畳み込み乗算 と、逆に下位から差し引かれる 負巡回畳み込み乗算 の他、$\alpha$ 倍されて巡回するパターンもある。1

通常 ”多倍長乗算に FFT 畳み込み乗算を利用する” という場合に用いている方法は正巡回畳み込み乗算であり、上位の演算結果が巡回しないようにあえて 0 を詰め込む技法が取られている。例えば

\begin{eqnarray}

\{x_i\} &=& \{3,1,4,1\}\\
\{y_i\} &=& \{5,9,2,6\}
\end{eqnarray}

の畳み込み乗算を考えると、

\begin{eqnarray}

&&(t^3+4t^2+t+3)(6t^3+2t^2+9t+5)\\
&=& 6t^6+26t^5+23t^4+61t^3+35t^2+32t+15
\end{eqnarray}

ということから、

\{3,1,4,1\}\times\{5,9,2,6\}=\{15,32,35,61,23,26,6\}

というのが畳み込み乗算結果になる。普通に DFT を経由してこの結果を出すには上位に $\{0\}$ を追加して 8 要素の DFT を行い、

\begin{eqnarray}

&&{\rm DFT}^{-1}({\rm DFT}(\{3,1,4,1,0,0,0,0\})*{\rm DFT}(\{5,9,2,6,0,0,0,0\}))\\
&=&\{15,32,35,61,23,26,6,0\}
\end{eqnarray}

ということになる…のだが、この演算を元のデータ数 4 要素のまま行うと

\begin{eqnarray}

&&{\rm DFT}^{-1}({\rm DFT}(\{3,1,4,1\})*{\rm DFT}(\{5,9,2,6\}))\\
&=&\{38,58,41,61\}
\end{eqnarray}

となる。この結果は多項式表現で言えば

\begin{eqnarray}

&&(t^3+4t^2+t+3)(6t^3+2t^2+9t+5) \pmod{t^4-1}\\
&=& 61t^3+(35+6)t^2+(32+26)t+(15+23)\\
&=& 61t^3+41t^2+58t+38
\end{eqnarray}

と等しい。これは 4 要素を超える部分が巡回して足し込まれた形になっている(多項式表現の 2 行目)。これを正巡回多項式という。逆に超えた部分が負になる

\begin{eqnarray}

&&(t^3+4t^2+t+3)(6t^3+2t^2+9t+5) \pmod{t^4+1}\\
&=& 61t^3+(35-6)t^2+(32-26)t+(15-23)\\
&=& 61t^3+29t^2+6t-8
\end{eqnarray}

の形を負巡回畳み込み乗算という。


離散荷重変換

さて、正巡回畳み込み乗算については通常の DFT を使って計算すればいいことは先述の通りなのだが、負巡回畳み込み乗算はどうするのかというと、離散荷重変換(Discrete Weighted Transform; DWT)なるものを使う。また新しい何かかというイメージもあるが、そんなに面倒くさいことでもない。

再度 DFT と正巡回畳み込み乗算の関係を復習しておくと

z_l = n^{-1}\sum_{k=0}^{n-1}

\left(\sum_{i=0}^{n-1}x_iw^{ik}\right)
\left(\sum_{j=0}^{n-1}y_jw^{jk}\right) w^{-kl}
= \sum_{i+j=l} x_iy_j + \sum_{i+j=l+n} x_iy_j

という感じである。実はこの式の最後の $\sum$ が巡回している部分を表してる。離散荷重変換ではここで新しい数列 $\{a_i\}$ を導入して

\begin{cases}

X_k = \sum_{j=0}^{n-1}a_j x_j w^{jk}\\
x_j = n^{-1}a_j^{-1} \sum_{k=0}^{n-1} X_k w^{-kj}
\end{cases}

という計算をする。入力となる数列 2 つ、出力となる数列を使って

\begin{eqnarray}

\{X_k\} &=& {\rm DWT}(\{x_j\}, \{a_j\})\\
\{x_j\} &=& {\rm DWT}^{-1}(\{X_k\}, \{a_j\})
\end{eqnarray}

と書こう。DFT との関係をまとめて書くなら

\begin{eqnarray}

\{X_k\} &=& {\rm DWT}(\{x_j\}, \{a_j\}) &=& {\rm DFT}(\{x_j\} * \{a_j\})\\
\{x_j\} &=& {\rm DWT}^{-1}(\{X_k\}, \{a_j\}) &=& \{a_j^{-1}\} * {\rm DFT}^{-1}(\{X_k\})
\end{eqnarray}

となるだろうか。

数式をモリモリ解いていくのは面倒だが、とりあえず $a_i=w^{i/2}$ とすれば

\begin{eqnarray}

\{Z_l\} &=& {\rm DWT}(\{x_i\}, \{a_i\}) * {\rm DWT}(\{y_j\}, \{a_j\})\\
\{z_l\} &=& {\rm DWT}^{-1}(\{Z_l\}, \{a_l\})\\
&=& \left\{\sum_{i+j=l} x_iy_j - \sum_{i+j=l+n} x_iy_j \right\}
\end{eqnarray}

つまり負巡回畳み込み乗算結果が求まることが分かる。もちろんこの原理は整数 DFT (NTT) でも使える。


効率化

元の記事で挙げた効率化と同様に元のデータや演算結果が実数であることを利用すると、複素数として無駄になっている部分にもデータを詰め込んで負巡回畳み込み乗算を効率的に求めることができる。元と同様に

\begin{eqnarray}

x'_j &=& x_{2j} + ix_{2j+1}\\
X_k &=& \sum_{j=0}^{n-1} x_j w^{j/2} w^{jk}\\
X'_k &=& \sum_{j=0}^{n/2-1} x'_j w^j w^{2jk}
\end{eqnarray}

とする($w$ は全体を通して 1 の原始 $n$ 乗根とする。また、$\{X'_k\}$ は項数 $n/2$ になっていることに注意。)と、 $\{X_k\}={\rm DWT}(\{x_j\},\{w^{j/2}\})$ と $\{X'_k\}={\rm DWT}(\{x'_j\},\{w^j\})$ の関係は

\begin{cases}

A'_k &= \dfrac{1}{2}(1+iw^{k+1/2})(X'_k - \overline{X'_{n/2-1-k}})\\
X_k &= X'_k - A'_k\\
X_{n/2-1-k} &= X'_{n/2-1-k} + \overline{A'_k}
\end{cases}

\begin{cases}

A_k &= \dfrac{1}{2}(1-iw^{-(k+1/2)})(X_k-\overline{X_{n/2-1-k}})\\
X'_k &= X_k - A_k\\
X'_{n/2-1-k} &= X_{n/2-1-k} + \overline{A_k}
\end{cases}

となる。よく見ると半分くらい計算されていないように見えるが、元データが実数であることを利用した共役の関係 ($X_{n-1-k}=\overline{X_k}$) で実際には求まっている。

これを一応コードっぽくすると以下のようになる。

void NegaRealDWT(Complex* X, int n, DFTDirection dir) {

// 注: 引数の n は元のデータ数の半分
const double t = M_PI / n;
if (dir == Backward) {
// Convert real to complex
for (int i = 0; i < n / 2; ++i) {
Complex a = (X[i] - X[n - 1 - i].conj()) * Complex(1 - sin(th * (i + 0.5)) , -cos(th * (i + 0.5))) / 2;
X[i] -= a;
X[n - 1 - i] += a.conj();
}
}
NegaComplexDFT(X, n, dir);
if (dir == Forward) {
// Convert complex to real
for (int i = 0; i < n / 2; ++i) {
Complex a = (X[i] - X[n - 1 - i].conj()) * Complex(1 - sin(th * (i + 0.5)) , cos(th * (i + 0.5))) / 2;
X[i] -= a;
X[n - 1 - i] += a.conj();
}
}
}

void NegaCyclicRealConvolution(double x[], double y[], int n, double z[]) {
Complex* X = static_cast<Complex*>(x);
Complex* Y = static_cast<Complex*>(y);
NegaRealDWT(X, n / 2, Forward); // n 個の実数 == n/2 個の複素数
NegaRealDWT(Y, n / 2, Forward);
for (int i = 0; i < n / 2; ++i)
Z[i] = X[i] * Y[i];
Complex* Z = static_cast<Complex*>(z);
NegaRealDWT(Z, n / 2, Backward);
}





  1. 多数桁計算における高速アルゴリズムの研究, 後保範, 早稲田大学大学院 博士論文