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

整数環FFT

適当な条件を満たす環や体では DFT を利用した時と同様の計算方法で 畳み込み乗算を行うことができる。 その際にはもちろん FFT と同様の式変換を行うことで $O(n\log n)$ の計算量で演算することが可能であり、 そのような方法を「整数環 FFT」もしくは 「数論変換 (Number-Theoretical Transform;NTT)」と呼ぶ。 特に、いくつか特殊な環や体で行う整数環 FFT については DGT (離散ガロア変換)、 DFT (離散フェルマー変換) といったように名前がつけられていることがある。

概要

適当な整数 $n$ に対して

  • 乗算における $n$ の逆元 $n^{-1}$ が存在する
  • 1 の原始 $n$ 乗根 $\omega$ が存在する

という条件を満たす環においては \[ A_k = \sum_{j=0}^{n-1} a_j \omega^{jk},\quad a_j = n^{-1} \sum_{k=0}^{n-1} A_k \omega^{-jk} \]

という相互変換が成り立つ。 これを用いて畳み込み乗算と同じ計算手順を踏むと同様に

\[ c_l = n^{-1}\sum_{k=0}^{n-1} \sum_{j=0}^{n-1} \sum_{j'=0}^{n-1} a_j b_{j'} \omega^{(j+j'-l)k} \] \[ \sum_{k=0}^{n-1}\omega^{jk} = \begin{cases} 0 & (j \bmod n \neq 0)\\ n & (j \bmod n = 0) \end{cases} \]

の 2 つが成り立つので、同様に畳み込み乗算を計算することができる。

64 bit 整数で行える剰余環 FFT

[JT03] では剰余環 $\mathbb{Z}/m\mathbb{Z}$ での計算やその見つけ方が紹介されている。 概要としてはまず $n$ を決め、$m=an+1$ として条件を満たすように 適当に決めた範囲で $a$ と $\omega$ を探索する。

ここでは $\log_2m$ が 25、51、64 付近になる組み合わせを探索した結果を 表1に示した。 $\log m\simeq 25, 51$ のものは倍精度浮動小数点数を用いて、 $\log m\simeq 64$ のものは 64 ビット整数を使って計算しやすい。 ただし $\log m\simeq 51, 64$ のものについては $a\cdot b \bmod m$ の計算を誤差なく行うために特殊な命令、もしくはコード上の工夫が必要となるので注意してほしい。

この表より小さな $n'$ で計算したい場合は倍数となる $n$ を用いて $w'=w^{n/n'}$ とすれば計算できる。

表1. 整数環FFTに使える数の例
$n$$\omega$$m$$\log_2m$
$2^{21}$38$11n+1$24.459
$3\times2^{20}$19$9n+1$24.755
$5\times2^{19}$19$10n+1$24.644
$2^{46}$67$27n+1$50.755
$3\times2^{45}$8$18n+1$50.755
$5\times2^{45}$13$7n+1$50.129
$2^{59}$87$27n+1$63.755
$3\times2^{58}$2$18n+1$63.755
$5\times2^{57}$51$19n+1$63.570

mod $2^{64}-2^{32}+1$

素数 $p=2^{64}-2^{32}+1$ を用いた剰余環 $\mathbb{Z}/p\mathbb{Z}$ は体となる。 そして 7 が原始 $p-1=2^{32}\cdot3\cdot5\cdot17\cdot257\cdot65537$ 乗根となるため $\omega=7^{(p-1)/n}$ を使うことで $n=2^{32}$ までの長さなら 2 基底 FFT ルーチンで対応できる。 さらに 1 回ずつまでであれば 3 基底や 5 基底ルーチンも適用させることでより長い $n$ でも計算することができる。

また、$7^{5(p-1)/192}=2$ が成り立つため、$n \lt 2^{26}$ に対して $\{a_j\}=\{(7^{5(p-1)/192/n})^j\}$ を利用することで整数演算を用いながら IBDWT と同じ原理を適用することも可能になる。

この体での計算では

\[ \begin{cases} 2^{64} = 2^{32} -1\\ 2^{96} = -1\\ 2^{128} = -2^{32}\\ 2^{192} = 1 \end{cases} \]

が成り立つので、$\bmod p$ の演算が容易になる。 具体的には一旦 128 bit になる計算結果を 32 bit 毎に区切って下の桁から順に $x_0$, $x_1$, $x_2$, $x_3$ とすると

\[ \begin{eqnarray} && x_0 + 2^{32}x_1 + 2^{64}x_2 + 2^{96}x_3\\ &=& 2^{64}x_2 + 2^{32}x_1 + (x_0-x_3) & (2^{96} = -1 \ {\rm を適用})\\ &=& 2^{32}(x_1+x_2) + x_0 - x_3 - x_2 & (2^{64} = 2^{32} - 1 \ {\rm を適用}) \end{eqnarray} \]

という形になり、演算命令として mod が必要なくなる。

Schönhage-Strassen Algorithm

FFT や NTT を利用した多倍長乗算の問題点の 1 つに、 1 変数に代入できる数の大きさに限界があることがあげられる。 一般的な倍精度浮動小数点数を使った計算では畳み込み演算結果が $2^{53}$ を超えると誤差が出る上に、 $2^{51}$ を超えるくらいでも整数化誤差が心配になる。 また、上記のような整数を使って計算誤差を無くした演算でも畳み込み演算結果が $m$ や $p$ を超えると計算を間違えてしまう。 Schönhage-Strassen Algorithm (SSA) では NTT を行う各要素も多倍長数にすることで そのようなオーバーフローを無くすことを提案している。

SSA は $n$ を元に $m=2^{n/2}+1$ を設定し $\mathbb{Z}/m\mathbb{Z}$ で計算する NTT である。 すると $\omega=2$ が原始 $n$ 乗根となるので上に挙げた条件を満たし計算ができる。 先に書いたとおり各要素を多倍長整数で実装する必要があるが、$\bmod{m}$ での実装なので $n/2$ bit ごとに区切り下位から交互に加算、減算していくことで計算ができる。

SSA のメリットとしては各要素が多倍長数になるのでワードサイズより小さくする分割が必要なく、 メモリ効率が高くなる点である。 また再帰的に多倍長演算を利用することができるので事実上計算桁数の上限が無いことも大きな利点である。 理論的な計算量は 2 段階の再帰計算で $O(n\log n\log^2n)$、 3 段階の再帰計算で $O(n\log n\log^2n\log^3n)$ と再帰回数が増える度に $\log^kn$ が増えていくが、実用的には 1 京桁程度の計算でも 2 段階の再帰計算で十分であることと $O(\log^3n)$ くらいからは実質的に定数と差がない程度に小さいため SSA の計算量は $O(n\log n\log^2n)$ とされている。