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

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

畳み込み乗算

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

通常 ”多倍長乗算に FFT を利用する” という場合に用いている 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}

についてもこれから述べる離散荷重変換を使うことで FFT と同様の計算量で計算ができます。

離散荷重変換

離散荷重変換(Discrete Weighted Transform; DWT)、また新しい何かかというイメージもありますが、FFT の前後に少し処理が入るだけの手法です。

再度 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}

という計算をします。この $\{x_j\}$, $\{a_j\}$, $\{X_k\}$ の関係を

\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}

となります。よく見ると半分くらい($n/2\lt k \lt n$)計算されていないように見えますが、元データが実数であることを利用した共役の関係 ($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. 多数桁計算における高速アルゴリズムの研究, 後保範, 早稲田大学大学院 博士論文 ↩