円周率.jp > 多数桁計算 > FFT > 整数環FFT

整数環FFT

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

Z/pZ の環において 自然数 ωωj \not\equiv 1 (0<jn), ωn≡1 が成り立つように pnω の組み合わせを選択すれば,複素数を用いる DFT とほぼ同じ式

Ak\sum_{j=0}%5e{n-1} aj ωjk,  ajn-1 \sum_{k=0}%5e{n-1} Ak ω-jk   (mod p)

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

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

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

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

ω'ω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