以下のソースはHSPコンテストに出した作品「FFTモジュール活用例２巨大整数の乗算」のソースです。
下の方にあるfftモジュールは著作権フリーなので自由に使って下さい
「fft」の使い方ですが、まずこの命令を使うには変数
「beki」
「LIST_SIZE」
「LIST_SIZE2」
「LIST_SIZE8」
の作成をしておいて下さい。
beki には高速フーリエ変換するサンプル数(2^n)のnを入れます。以下のソースで言うとサンプル数を8192個作りたかったので
beki=13
になっています。
LIST_SIZEには2^beki を入れて下さい
LIST_SIZE2はLIST_SIZE\2を、LIST_SIZE8にはLIST_SIZE×8を代入して下さい。
高速フーリエ変換したい実数(虚数)配列は当然ddim で作り、配列数はLIST_SIZEを指定して下さい
ここで注意ですがbeki には整数しか指定出来ませんので、変換する配列数は2のべき乗に限られます
命令「fft」の引数は３つで
・実数部(浮動小数点配列)
・虚数部(浮動小数点配列)
・逆変換か否か(1 or -1 -1で逆変換)
です。命令を実行すると実数部、虚数部に変換後の値が入ります。
下の方にあるfftモジュールは著作権フリーなので自由に使って下さい
「fft」の使い方ですが、まずこの命令を使うには変数
「beki」
「LIST_SIZE」
「LIST_SIZE2」
「LIST_SIZE8」
の作成をしておいて下さい。
beki には高速フーリエ変換するサンプル数(2^n)のnを入れます。以下のソースで言うとサンプル数を8192個作りたかったので
beki=13
になっています。
LIST_SIZEには2^beki を入れて下さい
LIST_SIZE2はLIST_SIZE\2を、LIST_SIZE8にはLIST_SIZE×8を代入して下さい。
高速フーリエ変換したい実数(虚数)配列は当然ddim で作り、配列数はLIST_SIZEを指定して下さい
ここで注意ですがbeki には整数しか指定出来ませんので、変換する配列数は2のべき乗に限られます
命令「fft」の引数は３つで
・実数部(浮動小数点配列)
・虚数部(浮動小数点配列)
・逆変換か否か(1 or -1 -1で逆変換)
です。命令を実行すると実数部、虚数部に変換後の値が入ります。
boxf
color 255,128,128
font "MS ゴシック",27,1
mes "Ａ"
color 255,255,255
pos 25,0
mes "× ＝"
color 128,255,128
pos 60,0
mes "Ｂ"
pos 140,2
button goto "答えを計算",*jikkou
color 255,128,128
boxf 0,30,314,480
color 128,255,128
boxf 314,30,640,480
pos 10,66
sdim k1,32768/2;(2^24)
mesbox k1,290,404,1,32766/2-1
sdim k2,32768/2;(2^24)
pos 332,66
mesbox k2,290,404,1,32766/2-1
ddim d1,32768/4
ddim d2,32768/4
ddim chasr,32768/4
ddim chasi,32768/4
stop
;計算実行
*jikkou
clrobj
color 255,255,255
boxf
beki=13;beki 13では8192配列が作成され1配列あたり4桁入るので最終出力が32768桁までなら計算できることになる
LIST_SIZE=1<<beki
LIST_SIZE2=LIST_SIZE>>1
LIST_SIZE8=LIST_SIZE<<3
strk1=strlen(k1)
strk2=strlen(k2)
hot4=4
repeat 32768/4/2;入力の最大桁は32768桁の半分
ccnt=LIST_SIZE2-1-cnt
sitenk1=strk1-cnt*4-4
if sitenk1<=0:hot4=strk1-cnt*4:sitenk1=0
d1.ccnt=double(strmid(k1,sitenk1,hot4));
if sitenk1<=0:break
loop
hot4=4
repeat 32768/4/2
ccnt=LIST_SIZE2-1-cnt
sitenk1=strk2-cnt*4-4
if sitenk1<=0:hot4=strk2-cnt*4:sitenk1=0
d2.ccnt=double(strmid(k2,sitenk1,hot4));
if sitenk1<=0:break
loop
fft d1,d2,1;フーリエ変換
memcpy chasr,d1,8*LIST_SIZE,0,0
memcpy chasi,d2,8*LIST_SIZE,0,0
ReTm=0.0
imTm=0.0
ReTN=0.0
imTN=0.0
ReTm=chasr
ImTm=chasi
ReTN=chasr
ImTN=chasi
d1=0.5*(Retm*ImTm+ReTN*ImTN)
d2=0.25*( ((ReTN+ImTN)*(ReTN-ImTN)) + ( (ImTm+ReTm)*(ImTm-ReTm) ) )
itjn=LIST_SIZE-1
repeat LIST_SIZE-1,1;フーリエ変換後の各要素同士を掛け算。この部分はこちらのページを参考に作った→http://homepage2.nifty.com/m_kamada/math/fftmul.htm
ReTm=chasr.cnt
ImTm=chasi.cnt
ReTN=chasr.itjn
ImTN=chasi.itjn
d1.cnt=0.5*(Retm*ImTm+ReTN*ImTN)
d2.cnt=0.25*( (ReTN+ImTN)*(ReTN-ImTN) + (ImTm+ReTm)*(ImTm-ReTm) )
itjn-
loop
ddim chasr,1
ddim chasi,1
fft d1,d2,-1;逆変換
sdim k2,4
sdim k1,LIST_SIZE*4+80
dim kekka,LIST_SIZE
w=LIST_SIZE-1;繰り上げ
w1=LIST_SIZE-2
kisuu=10000
repeat LIST_SIZE-1;繰り上げ
d1.w+=0.5
d1.w-=d1.w\1
if d1.w>=kisuu{
kuriaged=int(d1.w/kisuu)
d1.w-=kisuu*kuriaged
d1.w1+=kuriaged
}
kekka.w=int(d1.w+0.5)
w-
w1-
loop
d1.w+=0.5
d1.w-=d1.w\1
kekka.w=int(d1.w+0.5)
onf=0
w=0
w1=0
sdim k2,12;文字列に出力
repeat LIST_SIZE-1
if onf=0{
if kekka.cnt!0{
onf=1
k2=str(kekka.cnt)
w1=strlen(k2)
k2=strf("%0"+w1+"d",kekka.cnt)
}
}else{
w1=4
k2=strf("%04d",kekka.cnt)
}
memcpy k1,k2,w1,w,0
w+=w1
loop
color 0,0,0
pos 0,0:mes "答え＝"
mesbox k1,630,430,0,LIST_SIZE*4+80
wait 5
sdim kekka,1;メモリ解放
sdim k2,1
sdim d1,1
sdim d2,1
sdim chasr,1
sdim chasi,1
stop
#deffunc fft array A,array B,int inv
ddim c,LIST_SIZE
ddim d,LIST_SIZE
repeat beki;fft
bekii=LIST_SIZE>>(cnt+1)
repeat LIST_SIZE2
t2 = cnt\bekii;
t0 = (cnt/bekii)*bekii*2 + t2;
t1 = t0+bekii;
r4=A(t1)
i4=B(t1);
rm0=A(t0)-r4
im0=B(t0)-i4;
pii=3.14159265358979323846264338328*t2/bekii*inv;
dsin=sin(pii)
dcos=cos(pii);
A(t0) += r4;
B(t0) += i4;
A(t1) = rm0*dcos-im0*dsin;
B(t1) = rm0*dsin+im0*dcos;
loop
loop
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)
D(cnt)=B(ic)
loop
if inv=-1{
repeat LIST_SIZE
a.cnt=c.cnt/LIST_SIZE
b.cnt=d.cnt/LIST_SIZE
loop
}else{
memcpy a,c,LIST_SIZE8
memcpy b,d,LIST_SIZE8
}
return
