--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#const beki 3;2^bekiこの配列が作られます、例beki=23の時　約400万配列×400万配列が800万配列に格納されます
#const kisuu 10000;1～32768まで、浮動小数配列ひとつに格納できる上限数＝基数、何進数かという話、2^beki*log10底のkisuu が乗算後の値の桁数に値する

#module
#const beki beki@
#deffunc fft初期設定
n=1
repeat beki;桁数計算
n*=2
loop

ddim qcosq,n/2+1;sin cos予め計算格納
ddim qsinq,n/2+1
n2=n/2
repeat n/2+1
itjrad=3.1415926535897932384626*cnt/n2
qcosq.cnt=cos(itjrad)
qsinq.cnt=sin(itjrad)
loop
qcosq.n2=-1.0
qsinq.n2=0.0

itjn=1;バタフライインデックス作成
dim inind1,itjn
dim inind2,itjn*2
repeat beki
repeat itjn:inind2.cnt*=2:loop
dim inind1,itjn
memcpy inind1,inind2,4*itjn,0,0
dim inind2,itjn*2
memcpy inind2,inind1,4*itjn,0,0
repeat itjn:inind1.cnt+:loop
memcpy inind2,inind1,4*itjn,4*itjn,0
itjn*=2
loop
dim inind1,1;2の方に結果が
itjn=0

ddim chasr,n;キャッシュ
ddim chasi,n;キャッシュ
return

#deffunc fft array reitjo,array imitjo
q0q=1
q1q=n/2

repeat beki-1
qct0q=0
qct1q=q0q
qcn2q=0
repeat n/2
dup reitjoqccnt1,reitjo(inind2.qct0q)
dup imitjoqccnt1,imitjo(inind2.qct0q)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q
itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q
reitjoqccnt2=reitjoqccnt1-ita
imitjoqccnt2=imitjoqccnt1-itb
reitjoqccnt1+=ita
imitjoqccnt1+=itb
qct0q+
qct1q+
qcn2q+=q1q
if qct0q\q0q=0:qct0q+=q0q:qct1q+=q0q:qcn2q=0
loop
q0q*=2
q1q/=2
loop

qct0q=0
qct1q=q0q
qcn2q=0
repeat n/2
dup reitjoqccnt1,reitjo(inind2.cnt)
dup imitjoqccnt1,imitjo(inind2.cnt)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.cnt-imitjoqccnt2*qsinq.cnt
itb=imitjoqccnt2*qcosq.cnt+reitjoqccnt2*qsinq.cnt
chasr.cnt=reitjoqccnt1+ita
chasi.cnt=imitjoqccnt1+itb
chasr.qct1q=reitjoqccnt1-ita
chasi.qct1q=imitjoqccnt1-itb
qct1q+
loop

memcpy reitjo,chasr,8*n,0,0
memcpy imitjo,chasi,8*n,0,0
return



#deffunc ifft array reitjo,array imitjo
q0q=1
q1q=n/2

repeat beki-1
qct0q=0
qct1q=q0q
qcn2q=n2
repeat n/2
dup reitjoqccnt1,reitjo(inind2.qct0q)
dup imitjoqccnt1,imitjo(inind2.qct0q)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q
itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q
reitjoqccnt2=reitjoqccnt1+ita
imitjoqccnt2=imitjoqccnt1+itb
reitjoqccnt1-=ita
imitjoqccnt1-=itb
qct0q+
qct1q+
qcn2q-=q1q
if qct0q\q0q=0:qct0q+=q0q:qct1q+=q0q:qcn2q=n2
loop
q0q*=2
q1q/=2
loop

qct0q=0
qct1q=q0q
qcn2q=n2
repeat n/2
dup reitjoqccnt1,reitjo(inind2.cnt)
dup imitjoqccnt1,imitjo(inind2.cnt)
dup reitjoqccnt2,reitjo(inind2.qct1q)
dup imitjoqccnt2,imitjo(inind2.qct1q)
ita=reitjoqccnt2*qcosq.qcn2q-imitjoqccnt2*qsinq.qcn2q
itb=imitjoqccnt2*qcosq.qcn2q+reitjoqccnt2*qsinq.qcn2q
chasr.cnt=(reitjoqccnt1-ita)/n
chasi.cnt=(imitjoqccnt1-itb)/n
chasr.qct1q=(reitjoqccnt1+ita)/n
chasi.qct1q=(imitjoqccnt1+itb)/n
qct1q+
qcn2q-
loop


memcpy reitjo,chasr,8*n,0,0
memcpy imitjo,chasi,8*n,0,0
return

#global





screen 0,480,700
fft初期設定
randomize
n=1
repeat beki
n*=2
loop


ddim rex,n
ddim rey,n
ddim reitjo,n;実数代入先
ddim imitjo,n;虚数代入先
ddim chasr,n;キャッシュ
ddim chasi,n;キャッシュ

; gosub*indexsakusei
repeat n/2;「A＊B＝ｃ」のAとBにランダムな値を格納
rex.cnt=1.0*rnd(kisuu);2.7*sin(3.1415926535*490*cnt/n)+2.0*sin(3.1415926535*190.0*cnt/n)
loop
repeat n/2
rey.cnt=1.0*rnd(kisuu);1.0*sin(3.1415926535*50.0*cnt/n)+20.0*sin(3.1415926535*90.0*cnt/n)
loop

; aatim=gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)
; mes aatim
memcpy reitjo,rex,8*n,0,0
memcpy imitjo,rey,8*n,0,0

fft reitjo,imitjo

memcpy chasr,reitjo,8*n,0,0
memcpy chasi,imitjo,8*n,0,0


ReTm=0.0
imTm=0.0
ReTN=0.0
imTN=0.0
ReTm=chasr
ImTm=chasi
ReTN=chasr
ImTN=chasi
reitjo=(Retm*ImTm+ReTN*ImTN)/2.0
imitjo=( ((ReTN+ImTN)*(ReTN-ImTN)) + ( (ImTm+ReTm)*(ImTm-ReTm) ) )/4.0
itjn=n-1
repeat n-1,1
ReTm=chasr.cnt
ImTm=chasi.cnt
ReTN=chasr.itjn
ImTN=chasi.itjn
reitjo.cnt=0.5*(Retm*ImTm+ReTN*ImTN)
imitjo.cnt=0.25*( (ReTN+ImTN)*(ReTN-ImTN) + (ImTm+ReTm)*(ImTm-ReTm) )
itjn-
loop


ddim chasr,1;キャッシュ
ddim chasi,1;キャッシュ

ifft reitjo,imitjo


; mes gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)
; mes gettime(4)*3600000+gettime(5)*60000+gettime(6)*1000+gettime(7)-aatim

/*
repeat n
if absf(reitjo.cnt-int(reitjo.cnt))<0.7:if absf(reitjo.cnt-int(reitjo.cnt))>0.3:mes reitjo.cnt
; if absf(imitjo.cnt-int(imitjo.cnt))<0.7:if absf(imitjo.cnt-int(imitjo.cnt))>0.3:mes imitjo.cnt
loop
*/
;: bsave "reitjo",reitjo
; bsave "imitjo",imitjo
; stop
dim kekka,n;計算結果入るint型変数
gosub*繰り上げ計算
goto*hyouji2



*繰り上げ計算
w=n-1
w1=n-2
repeat n-1
kekka.w+=int(reitjo.w+0.5)
if kekka.w>=kisuu:kekka.w1+=kekka.w/kisuu:kekka.w\kisuu
w-
w1-
loop
kekka.w+=int(reitjo.w+0.5)
return

*hyouji;波形表示
color 0,0,0
repeat n
a=reitjo.cnt
if cnt=0:pos 0,int(240.0-a)
line cnt*640/n,int(240.0-a)
loop

color 255,0,0
repeat n
a=imitjo.cnt
if cnt=0:pos 0,int(240.0-a)
line cnt*640/n,int(240.0-a)
loop

color 0,255,0
repeat n
a=(imitjo.cnt*imitjo.cnt+reitjo.cnt*reitjo.cnt)
if cnt=0:pos 0,int(240.0-a)
line cnt*640/n,int(240.0-a)
loop
stop


*hyouji2
font "MS ゴシック",10,1
pos 0,0
mes "掛け算結果配列表示"
repeat limit(n,1,300)
if strlen(str(kekka.cnt))=1:mes "000"+str(kekka.cnt)
if strlen(str(kekka.cnt))=2:mes "00"+str(kekka.cnt)
if strlen(str(kekka.cnt))=3:mes "0"+str(kekka.cnt)
if strlen(str(kekka.cnt))>=4:mes str(kekka.cnt)
loop

pos 150,0
mes "掛ける数A配列表示"
repeat limit(n,1,300)
if strlen(str(int(rex.cnt)))=1:mes "000"+str(int(rex.cnt))
if strlen(str(int(rex.cnt)))=2:mes "00"+str(int(rex.cnt))
if strlen(str(int(rex.cnt)))=3:mes "0"+str(int(rex.cnt))
if strlen(str(int(rex.cnt)))>=4:mes str(int(rex.cnt))
loop

pos 300,0
mes "掛ける数B配列表示"
repeat limit(n,1,300)
if strlen(str(int(rey.cnt)))=1:mes "000"+str(int(rey.cnt))
if strlen(str(int(rey.cnt)))=2:mes "00"+str(int(rey.cnt))
if strlen(str(int(rey.cnt)))=3:mes "0"+str(int(rey.cnt))
if strlen(str(int(rey.cnt)))>=4:mes str(int(rey.cnt))
loop
stop

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
















これはhsp 3.31で動きます。他のverでの動作は保証できません

ソースの補足として、
命令として「fft」「ifft」「fft初期設定」の３つを登録した。それぞれ離散フーリエ変換、逆フーリエ、バタフライ演算のインデックス作成、となっている。
「fft」 「ifft」ともに、とる引数は実数配列（double）,虚数配列（double）で、変換結果がそれに返される。
配列の要素数は、2^beki個であることが条件、それ以上でも以下でも、１でもずれたらうまく計算されない。

実行すると、
実行画面の左には、A*B＝CのCの配列が表示されている
実行画面の真ん中には、A*B＝CのAの配列が表示されている
実行画面の右には、A*B＝CのBの配列が表示されている
bekiの値を増やすほど、多くの桁数の乗算が行える
仮にbekiを23にすれば、約400万配列×400万配列が800万配列に格納される
bekiが23でkisuuが10000の場合、約1600万桁×1600万桁=3200万桁の乗算となる


fftモジュールはもちろん乗算以外の用途にも使える。