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

整数環FFT

適当な条件を満たす環や体では DFT と同様の計算方法で畳み込み乗算を行うことができる。 その際にはもちろん FFT と同様の式変換を行うことで $O(n\log n)$ の計算量で演算することが可能であり、 そのような方法を「整数環 FFT」もしくは「NTT (Number-Theoretical Transform;数論変換)」と呼ぶ。 特に、いくつか特殊な環や体で行う整数環 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^{mk} = \begin{cases} 0 & (m \bmod n \neq 0)\\ n & (m \bmod n = 0) \end{cases} \]

も成り立つので、結局結果も畳み込み乗算になる。

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

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

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

表1. 整数環FFTに使える数の例
$n$$\omega$$m$$\log_2m$
$2^{17}$770$149n+1$24.2
$3\cdot2^{15}$26$173n+1$24.0
$5\cdot2^{14}$329$207n+1$24.0
$2^{35}$360$8319n+1$48.0
$3\cdot2^{33}$349$10857n+1$48.0
$3\cdot2^{32}$980$22396n+1$48.0
$5\cdot2^{31}$1942$25932n+1$48.0
$2^{50}$1521$1012n+1$60.0

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 の方法

FFT を利用した多倍長乗算の問題点の 1 つに、1 変数に代入できる数の大きさに限界があることがあげられる。 一般的な倍精度浮動小数点数を使った計算では畳み込み演算結果が $2^{53}$ を超えると誤差が出る上に、 $2^{51}$ を超えるくらいでも整数化誤差が心配になる。 また、上記のような整数を使って計算誤差を無くした演算でも畳み込み演算結果が $m$ や $p$ を超えると計算を間違えてしまう。

そこで、FFT を行う要素すら多倍長数にしてしまおうというのが Schönhage-Strassen の方法である。