なんかFFTよりも多倍長計算に向いている高速剰余変換（FMT）というものがあるらしい

で早速やってみた。

回転子であるωを4とし、すべての計算を257の剰余下で行う。
変換する全要素は8つであり、4^8を257で割った余りは1かつ、4^(8/2)を257で割った余りは-1(=256)であることがミソ

以下のソースでは4桁の10進数の数同士を乗算し結果を表示している



FFTを用いた計算よりもFMTが優れている点として、全て整数値同士の計算で完結するということがあげられる。
これは非常に重要だ。
FFTでは計算したい桁数を上げれば上げるほど、浮動小数点の誤差が蓄積してしまう。
doubleの精度を持ってしてでも1億桁×1億桁をやろうとすれば、double１要素につき4桁が限界だ。
これはdoubleつまり8byteに、10進数4桁≒1.66byteしか格納できなく、記憶するメモリ容量的にも無駄が多い。当然桁数あたりの必要な演算量も多くなってくる。
FMTでは4byte整数に約1.5～2byteは敷き詰められるしdoubleと違って演算スピードも早い。

誤差の話に戻るが、結果的にFFTを用いた多倍長の計算では1兆桁付近で精度の限界が生じてしまう。一方FMTは全て整数で計算するので理論上無限大まで誤差の蓄積なく計算することが可能だ。


というわけでFMT試作verのソースをUPする。
正直醜いスパゲッティなのであまり聞かないで・・・・
2014/2/3　YSRさんのご好意によりきれいなソースとなりました！ありがとうございます







#const beki 3
#const LIST_SIZE 1 << beki
#const LIST_SIZE2 LIST_SIZE >> 1
#const LIST_SIZE3 LIST_SIZE << 1
#const LIST_SIZE8 LIST_SIZE << 2
#const n LIST_SIZE
#const p 257
#const pn p/n
#const omega 4

#module
#deffunc fft array A, int inv
if(inv != 0){
dim c, LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =  cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic =  cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic = (iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a,c,LIST_SIZE8@
}
repeat beki@ ;fft
bekii = 1 << cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, t2*LIST_SIZE2@ / bekii, p@)
r4     = A(t1) * w \ p@
A(t1)  = (A(t0) + p@ - r4) \ p@
A(t0) += r4
A(t0)\ = p@
loop
loop
return

#deffunc ufft array A, int inv
repeat beki@ ;fft
bekii = LIST_SIZE2@ >> cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, t2 * LIST_SIZE2@ / bekii, p@)
r4     = A(t1)
A(t1)  = (A(t0) - r4 + p@) * w \ p@
A(t0) += r4
A(t0) \= p@
loop
loop
if(inv != 0){
dim c, LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic = cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic=(iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a, c, LIST_SIZE8@
}
return

#deffunc ifft array A,int inv
if(inv != 0){
dim c, LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic = cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic=(iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a, c, LIST_SIZE8@
}
repeat beki@ ;fft
bekii = 1 << cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, LIST_SIZE@ - t2 * LIST_SIZE2@ / bekii, p@)
r4     = A(t1) * w \ p@
A(t1)  = (A(t0) + p@ - r4) \ p@
A(t0) += r4
A(t0) \= p@
loop
loop
return

#deffunc iufft array A,int inv
repeat beki@;fft
bekii=LIST_SIZE2@>>cnt
repeat LIST_SIZE2@
t2     = cnt \ bekii
t0     = (cnt / bekii) * bekii * 2 + t2
t1     = t0 + bekii
w      = jouyo(omega@, LIST_SIZE@ - t2 * LIST_SIZE2@ / bekii, p@)
r4     = A(t1)
A(t1) =(A(t0) -r4 + p@) * w \ p@
A(t0) += r4
A(t0) \= p@
loop
loop
if(inv != 0){
dim c,LIST_SIZE@
repeat LIST_SIZE@ ;ビット逆順
ic  =cnt / 65536
ic  = (ic  & 0x00005555) << 1 | (ic  & 0x0000AAAA) >> 1
ic  = (ic  & 0x00003333) << 2 | (ic  & 0x0000CCCC) >> 2
ic  = (ic  & 0x00000F0F) << 4 | (ic  & 0x0000F0F0) >> 4
ic  = (ic  & 0x000000FF) << 8 | (ic  & 0x0000FF00) >> 8
iic = cnt \ 65536
iic = (iic & 0x00005555) << 1 | (iic & 0x0000AAAA) >> 1
iic = (iic & 0x00003333) << 2 | (iic & 0x0000CCCC) >> 2
iic = (iic & 0x00000F0F) << 4 | (iic & 0x0000F0F0) >> 4
iic = (iic & 0x000000FF) << 8 | (iic & 0x0000FF00) >> 8
ic=(iic * 32768 + ic / 2) >> (31 - beki@)
C(cnt) = A(ic)
loop
memcpy a, c, LIST_SIZE8@
}
return

#deffunc narasi array A,int inv
repeat LIST_SIZE@
arhgea  = a.cnt \ n@
a.cnt  /= n@
if(arhgea != 0) :a.cnt += (n@ - arhgea) * pn@ + 1
loop
return

#defcfunc jouyo int a, int inv, int cc
ninv=inv
m = 1.0 * cc
x = (1.0 * a) \ m
jouyoo = 1.0
repeat
/* n&1 はnが奇数(bitの1桁目が1)の時1、偶数(bitの1桁目が0)の時0を返す */
if(int(ninv) \ 2 == 1) :jouyoo = (jouyoo * x) \ m
ninv = ninv / 2
x = (x * x) \ m
if(ninv == 0) :break
loop
if(jouyoo < 0) :jouyoo += m
return int(jouyoo)
#global

dim aa, n
aa = 6,5,9,3,0,0,0,0
dim bb, n
bb = 7,8,6,5,0,0,0,0
dim cccc, n
dim dddd, n
memcpy cccc, aa, n*4
memcpy dddd, bb, n*4
;aa.1 *= -2
;aa.2 *= 4

ufft bb
ufft aa
repeat 8
bb.cnt *= aa.cnt
bb.cnt \= p
loop

ifft bb
narasi bb

repeat 7
bb.(cnt+1) += bb.cnt / 10
bb.cnt     \= 10
loop

pos 290, 0

repeat 4
mes cccc.cnt
pos ginfo_cx - 12, ginfo_cy - 18
loop

mes "*"
pos ginfo_cx-12, ginfo_cy-18

repeat 4
mes dddd.cnt
pos ginfo_cx - 12,ginfo_cy - 18
loop

mes "="
pos ginfo_cx - 12, ginfo_cy - 18

repeat 8
mes bb.cnt
pos ginfo_cx - 12, ginfo_cy - 18
loop

pos 0, 40
mes jouyo(400, 1024 * 64, 164 * 1024 * 128 + 1)
stop









ここで、FFTに時間間引きや周波数間引きの理論を取り入れることで、変換→逆変換の際2回あるビット逆順の計算を省くことが可能となる。
これも将来的にOpenCLへ落としこむ時、ランダムアクセスを減らしてパフォーマンス低下を防ぐ大きな意味を持つ。

さらにRADEON HD 7900シリーズ・・いわゆるGCNと言われるタイプのGPUコアは整数演算が非常に早い！
これでさらに面白くなってきそうだ！

もともとdouble精度の計算ができないGPUでも、mad24の機能で高速に整数演算ができる物が多いため、そういった点でもFFTよりFMTが有利だ！


というわけで次なるステップは、HSPCLへの移植だな


2019/10追記
FMTをPythonとCUDAで実装しました。ガウスルジャンドルで3億桁まで求めた記事をQiitaに公開しました。(この量だとlivedoorブログだと文字制限で圧倒的にオーバーする・・・)
https://qiita.com/Red_Black_GPGPU/items/e933e0d846b874b86c32
ソースはこちら
https://github.com/toropippi/FMT_QiitaSample