高速数論変換で多倍長整数の畳み込み乗算

次の素数 $p$¹ を法とする剰余環 $\mathbb{Z}/p\mathbb{Z}$ を考えます。

p = 2^{64}-2^{32}+1

性質

以下の自明な合同式を変形し、両辺に $2^{32}$ を次々に掛けると、

\begin{eqnarray*}
2^{64}-2^{32}+1 &\equiv& 0\mod{p} \\
2^{64} &\equiv& 2^{32}-1\mod{p} \\
2^{96} &\equiv& 2^{64}-2^{32}\mod{p} \\
&\equiv& 2^{32}-1-2^{32}\mod{p} \\
&\equiv& -1\mod{p} \\
2^{128} &\equiv& -2^{32}\mod{p} \\
2^{160} &\equiv& -2^{64}\mod{p} \\
&\equiv& -2^{32}+1\mod{p} \\
2^{192} &\equiv& -2^{64}+2^{32}\mod{p} \\
&\equiv& -2^{32}+1+2^{32}\mod{p} \\
&\equiv& 1\mod{p} \\
\end{eqnarray*}

たとえば4つの32ビット整数 $x_0,x_1,x_2,x_3$ からなる128ビット整数 $X$ を $p$ で割った余りは、

\begin{eqnarray*}
X &=& x_0+x_12^{32}+x_22^{64}+x_32^{96} \\
&\equiv& x_0+x_12^{32}+x_2\left(2^{32}-1\right)+x_3\left(-1\right)\mod{p} \\
&\equiv& x_0-x_2-x_3+\left(x_1+x_2\right)2^{32}\mod{p} \\
\end{eqnarray*}

のように変形することで、剰余算なしで計算することが出来ます。
次節の数論変換で剰余算を多用するため、この性質からは大きな恩恵を受けられます。

数論変換

$p$ の原始 $p-1$ 乗根は $7$ なので、フェルマーの小定理より以下が成り立ちます。

7^{p-1}\equiv{1}\mod{p}

$p-1$ を因数分解すると、

\begin{eqnarray*}
p-1 &=& 2^{64}-2^{32} \\
&=& 2^{32}\left(2^{32}-1\right) \\
\end{eqnarray*}

以下のような $r^{[c]}_n$² を定義すると、

\left(\left(7^{2^{32}-1}\right)^{2^{\left(32-c\right)}}\right)^{n+1}\equiv{r^{[c]}_n}\mod{p}\qquad{%
\forall{\left(c,n\right)}\in\mathbb{Z},0<c\leq{32},0\leq{n}<2^c
}

次の等式が成り立ちます。

\left\{
\begin{eqnarray*}
r^{[c]}_a r^{[c]}_b &\equiv& r^{[c]}_{a+b+1}\mod{p} \\
r^{[c]}_{2^ca+b} &\equiv& r^{[c]}_b\mod{p} \\
r^{[c]}_{2^{c-1}a+b}+r^{[c]}_b &\equiv& 0\mod{p} \\
r^{[c]}_{a-2} r^{[c]}_{2^c-a} &\equiv& 1\mod{p} \\
\end{eqnarray*}
\right.
\quad
\forall{\left(a,b\right)}\in\mathbb{Z}

この $r^{[c]}_n$ を用いて、入力 $f(x),g(x)$ に対して下記のように畳み込み乗算 $(f*g)$ を計算することができます。

\left\{
\begin{eqnarray*}
F^{[c]}_j &=& \sum_{i=0}^{2^c-1}{%
r^{[c]}_{ij}f(i)
} \\
G^{[c]}_j &=& \sum_{i=0}^{2^c-1}{%
r^{[c]}_{ij}g(i)
} \\
\left(F*G\right)^{[c]}_j &=& F^{[c]}_j\times{G^{[c]}_j} \\
\left(f*g\right)^{[c]}_j &=& \mathrm{inv}_p\left(2^c\right)\sum_{i=0}^{2^c-1}{%
\mathrm{inv}_p\left(r^{[c]}_{ij}\right)\left(F*G\right)^{[c]}_j
} \\
\end{eqnarray*}
\right.

なお $\mathrm{inv}_p(x)$ は $x\times\mathrm{inv}_p(x)\equiv{1}\mod{p}$ を満たす値、つまりは $\mathbb{Z}/p\mathbb{Z}$ における逆元です。

$r^{[c]}_n$ の $\mathbb{Z}/p\mathbb{Z}$ における対称性を利用して高速フーリエ変換と同様のバタフライ演算³をおこなうことで乗算回数を $\mathcal{O}\left(n\log{n}\right)$⁴ に減らすことができます。

さらに乗算回数を減らす

ここまでは一般的な高速数論変換の話でした。
ここからは $r^{[6]}_n$ の世界での話になります。

なぜ $r^{[6]}_n$ かというと $7^{2^{26}\left(2^{32}-1\right)}\equiv{r^{[6]}_0}\equiv{2^{39}}\mod{p}$ なので、任意の $n\in\mathbb{Z}$ について $r^{[6]}_n$ を2の冪乗で表すことができるからです。
これはすなわち $r^{[6]}_n$ の乗算をビットシフト⁵に置き換えられることを意味します。

実装例

最後に実装例⁶へのリンクを置いておきます。

1. 10進数で $\texttt{18446744069414584321}$
   16進数で $\texttt{0xffffffff00000001}$ ↩

2. 本記事では、右肩の角括弧で囲った ${^{[c]}}$ の $c$ は、冪乗ではなく要素数の対数となるような変数を表しています。 ↩

3. 実装例では再帰せずにループで処理しています。 ↩

4. 実装例では $n\left(2+3\log_2{n}\right)$ 回の乗算に相当する計算のうち $r^{[6]}_n$ の乗算部分をビットシフトに置き換えたことで $2n$ 回まで減っています。 ↩

5. 実際には $p$ での剰余を求めなければならないので加算・減算・比較・ビットシフト演算が必要です。 ↩

6. 本記事執筆時点で、繰り上がりの計算と64bit⇔63bitの変換の実装が正しくないため、出力結果の積は正しくありません。 計算結果は正しくなりましたが、どこかで桁溢れが発生するようで、おおよそ $29\mathrm{bits}\times14\texttt{~}15\texttt{桁}$ 程度の多倍長整数同士までしか正しく計算できません。バタフライ演算の実装を手抜きしたのが良くなかったのかな… ↩