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

整数環FFT

整数環FFTは,剰余環 $\mathbb{Z}/p\mathbb{Z}$ でFFTを行う方式で, NTT(Number-Theoretical Transform;数論変換) も呼ばれる. 特に法に取る値 $p$ により DGT(Discrete Galois Transform;離散ガロア体変換), FFT(Fast Fermat Transform;高速フェルマー変換) といった特別な名がある.

$\mathbb{Z}/p\mathbb{Z}$ の環において自然数 $\omega$ が $\omega^j \neq 1$ (for $0\lt j\lt n$),$\omega^n=1$ を満たすように $p$, $n$, $\omega$ の組み合わせを選択すれば, 複素数を用いる DFT と(ほぼ)同じ式

\[ A_k \equiv \sum_{j=0}^{n-1} a_j \omega^{jk},\quad a_j \equiv n^{-1} \sum_{k=0}^{n-1} A_k \omega^{-jk} \pmod{p} \]

が成り立つ. この式から,DFT から FFT を導出したものと同じ変形ができるので, $O(n\log n)$ の計算量で変換可能なこと,畳込み乗算を実現できること, などが示せる.

簡単なところでは $p=\omega^{n/2} + 1$ と定義すれば条件を満たすが,$n$ が少し大きくなると $\omega=2$ でもすぐに $p$ が 1 変数で表せない数になるので,この場合は 1 つ 1 つの要素 ($a_i$) を多倍長数で実装する必要がある.

そのため,計算しようとする $n$ と $p$ の大きさに合わせて適当な $\omega$ を事前に計算しておくとよい. [JT03] にある値にいくつか追加し, 表1に示してみた. ここに示した値は $\log p \simeq 24, 48$ のものは倍精度浮動小数点数を用いて, $\log p \simeq 60$ のものは 64 ビット整数を使って計算できる例である. $\log p \simeq 48, 60$ のものについては $a b\bmod p$ の計算を誤差なく行わなければならないため, 特殊な命令,もしくはコード上の工夫が必要となる.

なお,表 1 にある $n$ より小さな $n'$ (ただし $n' | n$) での計算を行いたい場合,回転子 $\omega'$ を

\[ \omega' = \omega^{n / n'} \]

と再定義すれば良い.

表1. 整数環FFTに使える数の例
n ω p log2p
217 770 149n+1 24.2
3・215 26 173n+1 24.0
5・214 329 207n+1 24.0
235 360 8319n+1 48.0
3・233 349 10857n+1 48.0
3・232 980 22396n+1 48.0
5・231 1942 25932n+1 48.0
250 1521 1012n+1 60.0