FMT(高速剰余変換)
2002/12/20 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
2001/6/9 作成のものをweb化
1. 概要
FMT(Fast Modulo Transform, 高速剰余変換)はFFT(高速フーリエ変換)の原理に剰余(mod)
の拡張を追加して作成したものである。FFTとFMTの計算式は同一であるが、与える係数の
順が異なる。FMTは順変換の結果が剰余として意味つけられて多数桁計算への応用ではFFT
より見通しが良い。このため、連立一次方程式や対称行列の固有値問題の厳密解の高速計算
への応用が期待できる。
FMTの詳細はFMT原理 | 、
FMT応用 | 及び
FMTを利用した多数桁計算プログラム | を参照。
FMTを利用した円周率の値の計算はπ計算プログラム | を参照。
2. FMTの原理
中国剰余定理(百五減算)で使用するn個の素数Pkの代りにPk
=E-ωkを利用する。
ここで、Eは多項式を形成する記号で,ωは1の原始n乗根である。
更に、mod(f,E-ωk)をEの多項式関数fの記号Eを数値ωkで
置きかえるものと定義する。
このとき、FMT順変換とFMT逆変換を下記のように定義する。
(1) FMT順変換
f=an-1En-1+...+a1E+a0からfk
=mod(f,E-ωk)を求める。
fk=an-1ω(n-1)k+...+a1ωk
+a0 , k=0,1,...,n-1
(2) FMT逆変換
gk=mod(g,E-ωk)からもとの
g=bn-1En-1+...+b1E+b0の係数
bkを求める。
bk=(gn-1ω-(n-1)k+...+g1ω-k
+g0)/n , k=0,1,...,n-1
3. FMTの計算方法
(1) FMTによる畳込み演算
Eの多項式関数fとgの畳込み演算h=f・gをFMTで計算する。ここで、
f=an-1En-1+...+a1E+a0,
g=bn-1En-1+...+b1E+b0で
h=c2n-1E2n-1+c2n-2E2n-2+...+c1E
+c0とし、ωは1の原始2n乗根とする。
(a) 順変換 :fk=mod(f,E-ωk)=an-1ω(n-1)k+
...+a1ωk+a0 , gk=mod(g,E-ωk)
=bn-1ω(n-1)k+...+a1ωk+b0
(b) 項別乗算 :hk=fk・gk
(c) 逆変換 :hkから畳込み演算結果h=f・gの係数ckを計算。ここで、
k=0,1,...,2n-1
ck=(h2n-1ω-(2n-1)k+...+h1ω-k
+h0)/2n
(2) FMTとFFTの違い
FMTとFFTは計算式は同一であるが、下記が異なる。
(a) FMTは高位要素(上位桁)から順に与える、FFTは下位要素(下位桁)から順に与える。
多項式係数ak,k=0,1,...,n-1と計算用配列A(n)の対応は下記のようになる。
FFT: A(k+1)=ak , FMT: A(k+1)=an-k
(b) 多数桁乗算でFMTは桁合わせが一致、FFTでは下位桁がゼロとなり、1要素分シフト要。
(c) FFTのωは複素数(ω=e2πi/n)に限定されるが、FMTは整数にも拡張。
(3) mod(En-1)上でのFMTの計算方法
nが2のベキ乗の場合の例で示す。ωのテーブルTB(k)=ωn-k-1とする。
逆変換はTB(k)=ω-(n-k-1)
下記の計算を(m,l)=(n/2,1),(n/4,2),...,(1,n/2)の組に対して行い、各組ごとにWをAに入れ戻す。
dimension A(m,2,l),W(m,l,2),TB(m,l)
do j=1,l
do k=1,m
WA = A(k,1,j)*TB(m,j)
W(k,j,1) = A(k,2,j) - WA
W(k,j,2) = A(k,2,j) + WA
end do
end do
(4) mod(En-ωn/p)上でのFMTの計算方法
mod(En-1)上でのFMTの計算に下記を追加する。
(a) 順FMT:入力係数akにωk/pを乗算し、ak=ak
・ωk/pと変換する。
(b) 逆FMT:出力係数hkにω-k/pを乗算し、hk=hk
・ω-k/pと変換する。
ここで,k=0,1,...,n-1, p=0,1,...,n/2-1