ab.com コミュニティ

by **通りすがり** » 2007年7月02日(月) 21:45

ソースをよく読んでないのですが、FFT前の実数の後半半分はゼロで埋める必要があります。

http://www.nextftp.com/swlabo/m0_pctech ... _815_1.htm

by **のぶあや** » 2007年2月06日(火) 21:39

多倍長整数同士の乗算をFFT(高速フーリエ変換)を用いてやろうとしているのですが、正しく計算させることができません。
「多倍長整数 FFT」などでググってコーディングしたのですが、もしかしたら手順が違っているのかもしれません。

多倍長整数の構造

[ここをクリックすると内容が表示されます]

乗算の手順

[ここをクリックすると内容が表示されます]

多倍長整数の乗算

[ここをクリックすると内容が表示されます]

コード: 全て選択

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

FFTクラス(C言語による最新アルゴリズム事典から一部変更して移植)

[ここをクリックすると内容が表示されます]

参考元：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

ab.com コミュニティ

多倍長整数のFFT乗算について

返信する

展開ビュートピックのレビュー: 多倍長整数のFFT乗算について

実数の後半にゼロ埋めが抜けているのでは？

多倍長整数のFFT乗算について

多倍長整数のFFT乗算について

返信する

展開ビュー トピックのレビュー: 多倍長整数のFFT乗算について

実数の後半にゼロ埋めが抜けているのでは？

多倍長整数のFFT乗算について

展開ビュートピックのレビュー: 多倍長整数のFFT乗算について