「多倍長整数 FFT」などでググってコーディングしたのですが、もしかしたら手順が違っているのかもしれません。
多倍長整数の構造
[ここをクリックすると内容が表示されます]
乗算の手順・変数を「Dim SN As *DWord」として「calloc」でメモリを確保
・SN[0]： abbb bbbb bbbb bbbb bbbb bbbb bbbb bbbb （2進表記）
a:正数なら「0」、負数なら「1」
b:配列の要素数を記録（SN=calloc(SizeOf(DWord)*6)なら「5」）
・数値は10000進数で小さい桁から記録
「123456789」という数値ならSN[1]=6789, SN[2]=2345, SN[3]=1
・SN[0]： abbb bbbb bbbb bbbb bbbb bbbb bbbb bbbb （2進表記）
a:正数なら「0」、負数なら「1」
b:配列の要素数を記録（SN=calloc(SizeOf(DWord)*6)なら「5」）
・数値は10000進数で小さい桁から記録
「123456789」という数値ならSN[1]=6789, SN[2]=2345, SN[3]=1
[ここをクリックすると内容が表示されます]
多倍長整数の乗算
乗算する2数の大きいほうが10000進数でN桁の場合
・配列の要素を2N以上の２のべきにする
・2数をそれぞれフーリエ変換する
・実部の同じ位置にあるものを乗算する
・乗算した配列を逆フーリエ変換する
・10000進数になるように桁上がり処理をする
・配列の要素を2N以上の２のべきにする
・2数をそれぞれフーリエ変換する
・実部の同じ位置にあるものを乗算する
・乗算した配列を逆フーリエ変換する
・10000進数になるように桁上がり処理をする
[ここをクリックすると内容が表示されます]
FFTクラス(C言語による最新アルゴリズム事典から一部変更して移植)
コード: 全て選択
TypeDef SNUM = *DWord
Dim SN As SuperNumber '多倍長整数クラス
Dim SN_A As SNUM, SN_B As SNUM
'多倍長整数を作成
SN_A=SN.CreateSN(3)
SN_B=SN.CreateSN(3)
'数値を入力（文字列で入力する関数があるが分かりやすくするために省略）
'SN_A= 111,111,111
SN_A[1]=1111 '1～4桁
SN_A[2]=1111 '5～8桁
SN_A[3]=1 '9～12桁
'SN_B= 111,111
SN_B[1]=1111 '1～4桁
SN_B[2]=11 '5～8桁
'数値の確認
MessageBox(0, SN.Str(SN_A), "", MB_OK)
MessageBox(0, SN.Str(SN_B), "", MB_OK)
'乗算
SN.by(SN_A, SN_B)
'演算結果の確認
MessageBox(0, SN.Str(SN_A), "", MB_OK)
'終了
End
Class SuperNumber
Public
/*****************
数値生成・破棄
*****************/
Function CreateSN(figure As DWord) As SNUM
Dim buf As *DWord
buf=calloc(4*(figure+1))
buf[0]=figure
CreateSN=buf
End Function
/*****************
数値出力・表現
*****************/
'クラス内変数に文字列を格納し、ポインタで返す
Private
_Str_buf As BytePtr'コンストラクタで0初期化
Public
Function Str(SN As SNUM) As BytePtr
free(_Str_buf)
If SN[0]=0 then _Str_buf=0: Str=_Str_buf: Exit Function
Dim figure As DWord, c As Long, ctr As Long
figure=GetFigure(SN)
_Str_buf=calloc(figure+2)
ctr=figure\4: If (figure Mod 4)<>0 then ctr=ctr+1
If (SN[0] AND &h80000000)=&h80000000 then _Str_buf[0]=Asc("-")
wsprintf(_Str_buf+lstrlen(_Str_buf), "%1u", SN[ctr])
For c=ctr-1 to 1 Step -1
wsprintf(_Str_buf+lstrlen(_Str_buf), "%04u", SN[c])
Next
Str=_Str_buf
End Function
'格納されている桁数(メモリサイズではない)を求める
Function GetFigure(SN As SNUM) As DWord
If SN[0]=0 then GetFigure=0:Exit Function
Dim c As Long, ctr As DWord
ctr=0
For c=1 to (SN[0] AND &h7FFFFFFF)
If SN[c]<>0 then ctr=c
Next
If ctr=0 then GetFigure=0 :Exit Function
If SN[ctr]>=1000 then: c=4
Elseif SN[ctr]>=100 then: c=3
Elseif SN[ctr]>=10 then: c=2
Else: c=1
End If
GetFigure=4*(ctr-1)+c: Exit Function
End Function
/***********
演算処理
***********/
cl As FFT_CooleyTukey 'FFTクラス
Sub by(ByRef SN_A As SNUM, SN_B As SNUM)
Dim len As DWord, c As DWord
'len*2以上の2のべき
If SN_A[0]=>SN_B[0] then len=SN_A[0] else len=SN_B[0]
c=1: len=len*2
Do
c=c*2
If c>=len then len=c:Exit Do
Loop
'Double型に変換
Dim ar As *Double, ai As *Double, br As *Double, bi As *Double, ci As *Double
ar=calloc(SizeOf(Double)*len):ai=calloc(SizeOf(Double)*len)
br=calloc(SizeOf(Double)*len):bi=calloc(SizeOf(Double)*len):ci=calloc(SizeOf(Double)*len)
For c=1 to SN_A[0]
ar[c-1]=SN_A[c] As Double
Next
For c=1 to SN_B[0]
br[c-1]=SN_B[c] As Double
Next
'フーリエ変換
cl.FFT(len, ar, ai, FALSE)
cl.FFT(len, br, bi, FALSE)
'同じ位置同士を乗算
For c=0 to len-1
ar[c]=ar[c]*br[c]
Next
'逆フーリエ変換
cl.FFT(len, ar, ci, TRUE)
'桁上がり処理
Dim _SN As SNUM
Dim carry As DWord, buf As Double
_SN=CreateSN(len) '新規に作成
carry=0
For c=1 to len
buf=ar[c-1]+carry
_SN[c]=(buf Mod 10000)
carry=buf\10000
Next
free(SN_A): SN_A=_SN '新規に作成したものに交換
free(ar):free(ai):free(br):free(bi):free(ci)
End Sub
End Class
[ここをクリックすると内容が表示されます]
参考元：http://oku.edu.mie-u.ac.jp/~okumura/algo/
コード: 全て選択
Class FFT_CooleyTukey
Private
'FFT関数内のStatic変数
Last_n As Long
BitRev As *Long
SinTbl As *Double
/*** 関数{\tt fft()}の下請けとして三角関数表を作る. ***/
Sub Make_SinTbl(n As Long, SinTbl As *Double)
Dim i As Long, n2 As Long, n4 As Long, n8 As Long
Dim c As Double, s As Double, dc As Double, ds As Double, t As Double
n2=n/2: n4=n/4: n8=n/8
t=Sin(_System_PI/n)
dc=2*t*t: ds=Sqr(dc*(2-dc))
t=2*dc
c=1: SinTbl[n4]=1
s=0: SinTbl[0]=0
i=1: While i<n8
c=c-dc: dc=dc+t*c
s=s+ds: ds=ds-t*s
SinTbl=s: SinTbl[n4-i]=c
i=i+1: Wend
If n8<>0 then SinTbl[n8]=Sqr(0.5)
i=0: While i<n4: SinTbl[n2-i]=SinTbl: i=i+1: Wend
i=0: While i<n2+n4: SinTbl[i+n2]=(-SinTbl): i=i+1: Wend
End Sub
/*** 関数{\tt fft()}の下請けとしてビット反転表を作る. ***/
Sub Make_BitRev(n As Long, BitRev As *Long)
Dim i As Long, j As Long, k As Long, n2 As Long
n2=n/2: i=0: j=0
Do
BitRev=j
i=i+1
If i>=n then Exit Do
k=n2
While k<=j: j=j-k: k=k/2: Wend
j=j+k
Loop
End Sub
Public
Function FFT(n As Long, x As *Double, y As *Double, inv As Long)
Dim i As Long, j As Long, k As Long, ik As Long, h As Long
Dim d As Long, k2 As Long, n4 As Long', inv As Long
Dim t As Double, s As Double, c As Double, dx As Double, dy As Double
/* 準備 */
n4=n/4
If ((n<>Last_n) or (n=0)) then
Last_n=n
free(SinTbl):SinTbl=0
free(BitRev):BitRev=0
If n=0 then FFT=0:Exit Function
SinTbl=calloc((n+n4)*SizeOf(Double))
BitRev=calloc(n*SizeOf(Long))
If ((SinTbl=0) or (BitRev=0)) then FFT=1: Exit Function
Make_SinTbl(n, SinTbl)
Make_BitRev(n, BitRev)
End If
i=0: While i<n
j=BitRev
If i<j then
t=x: x=x[j]: x[j]=t
t=y: y=y[j]: y[j]=t
End If
i=i+1
Wend
k=1: While k<n
h=0: k2=k+k: d=n/k2
j=0: While j<k
c=SinTbl[h+n4]
If inv=TRUE then s= (-SinTbl[h]) else s=SinTbl[h]
i=j: While i<n
ik=i+k
dx=s*y[ik]+c*x[ik]
dy=c*y[ik]-s*x[ik]
x[ik]=x-dx: x[i]=x[i]+dx
y[ik]=y[i]-dy: y[i]=y[i]+dy
i=i+k2: Wend
h=h+d
j=j+1: Wend
k=k2: Wend
If inv<>TRUE then
i=0:While i<n
x[i]=x[i]/n: y[i]=y[i]/n
i=i+1:Wend
End If
End Function
End Class