Assembler de FFT

99/5/14

FPU は恐い

最近のコンパイラは賢いためかなり最適化されたコードを吐いてくれますが、それでもやはりアセンブラを使った方が(よっぽど間抜けなコードを書かない限り)一般には速くなります。しかし幾つか例外があります。その一つは FPU 命令を使った場合で、要するに浮動少数演算を扱う場合です。

詳しくはここを読んで頂くとして、具体的に FPU の何が問題かと言うと、主に

FPU レジスタはスタック形式であり、常に現在の状態を把握する必要がある。
ペアリングの制限が厳しい。
命令の制御が戻るタイミングと計算結果が使用できるようになるタイミングがずれている。
intel 系と AMD 系の FPU 命令の性質が異なるため、どちらか一方に合わせた最適化をすると大抵はもう片方が遅くなる。

という事があげられます。FPU 命令にはこれらの制約があるために手動で最適化しようとすると、1 つや 2 つの FPU 命令ならともかく、何十個も FPU 命令が続くような場合はかなり難しいです。てゆーか死にます。それで結局浮動少数を扱う場合には、少々の無駄があったとしても、下手にアセンブラで書くより C で書いてコンパイラに最適化してもらった方が速い場合が多いのです。

もちろん玉砕覚悟で FPU 命令の手動最適化をすればアセンブラの方が速くなりますが、暗号のクラックコンテストの様なコンテスト用のプログラムでもない限りあまりやる気がしませんね。FPU の最適化に費やす時間があったら他の部分のコードを書いていた方が良いです。

ここで具体例(失敗例 (^^;)をあげてみましょう。以前作った高速フーリエ変換のプログラムがどうも遅い様な気がしたので演算部分をアセンブラで書き直したのですが、面倒くさかったので FPU の最適化をしないで馬鹿正直なコードを書いたら逆に遅くなりました(笑)。手動で書いたコードよりもコンパイラが吐いたコードの方が 1.2 倍程速かったです。うお、コンパイラに負けた(^^;

下におまけで載せている VC++6.0 用のコードを見ていただくと分かるように、フーリエ変換の様な計算では FPU 命令が山の様に出てきます。こんなものを手動で最適化する気は起きません。誰か最適化して下さい(笑い)。という訳で、今度から FPU 命令を使うような場合はアセンブラ化しないことにしました(駄目杉)。ちなみに、固定少数を使うならばアセンブラ化した方がはるかに速くなります( 精度を気にしなければ MMX で同時演算も行えますし)。

(おまけ: パワースペクトル計算ルーチン)

重みの sin 値と cos 値(double 型配列、サイズ (int)(log(データ数)/log(2)+1)*データ数/2) とビット反転用の最大ビット数を PreCalcPowerSpec() で計算してから、 data に COMPLEX 型のデータを入力し(必ず 2 の倍数個)、データと重みを CalcPowerSpec() に渡すとパワースペクトルが data の実数パートに(半分だけ)入って戻る。

// 複素型構造体   
typedef struct
{
 double r;  // 実数パート
 double i;  // 虚数パート
}COMPLEX;


//------------------------------------
// パワースペクトル計算の前準備
//
// バタフライ計算の sin と cos の重みテーブル計算	   
// ビットの反転で使用する最大 BIT の計算をする
BOOL PreCalcPowerSpec(double* lpdCos,  // cos テーブル
					double* lpdSin,  // sin テーブル
					LPDWORD lpdwMaxBit, // 最大ビット数
					DWORD dwFftPoint  // ポイント数
					){

	DWORD i,i2;
	double dRad;
	DWORD dwStage = (DWORD)(log(dwFftPoint)/log(2));

	// Sin , Cos テーブルセット
	for(i=1;i< =dwStage;i++){
		for(i2=0;i2< dwFftPoint/2;i2++){ 
		dRad = 2.*PI*(double)(i2)/(double)(1<< i); 
		lpdCos[i*dwFftPoint/2+i2] = cos(dRad);				 
		lpdSin[i*dwFftPoint/2+i2] = sin(dRad);		
		}	 
	}

	// 最大 BIT の計算 
	i=dwFftPoint-1; 
	*lpdwMaxBit = 0;  
	while(i>0){i>>=1;(*lpdwMaxBit)++;};

	return TRUE;
}




//------------------------------------
// パワースペクトル計算のアセンブラ版
// FPU に手を出すとろくな目に会わない(--;
VOID CalcPowerSpec(COMPLEX* data,// FFT を行うデータ
				   DWORD dwFftPoint, // FFT の点数
				   DWORD dwMaxBit,  // dwFftPoint のビット数
				   double* dCosTable,// Cos のテーブル
				   double* dSinTable// Sin のテーブル
				   ){

	// FFT 用変数
	DWORD dwStage;
	DWORD dwPart;
	DWORD dwDataHeight,dwDataHeight2;

	// データのポインタ
	DWORD dwDataPos;

	// 雑用
	DWORD i,j,k;
		
	//------------
	// 計算部

	 dwPart = 1;
	 dwStage = (DWORD)(log(dwFftPoint)/log(2));

	 _asm{

		 //  for(dwStage = (DWORD)(log(dwFftPoint)/log(2));dwStage>0;dwStage--)
L_STAGE:

	 // dwDataHeight = (DWORD)1 <<  dwStage; // 2^dwStage
	 // dwDataHeight2 = dwDataHeight/2;

	 mov eax,1
	 mov ecx,dwStage
	 shl eax,cl
	 mov dwDataHeight,eax

	 mov ebx,eax
	 shr ebx,1
	 mov dwDataHeight2,ebx

      //--------------------------------
      // バタフライ演算開始

	  //    for(j = 0; j <  dwPart ;j++){
	  xor edi,edi // edi = j

L_PART:

	 // for(k = 0 ; k <   dwDataHeight2 ; k++){
	 xor esi,esi // esi = k

	 // dwDataPos1 = dwDataHeight*j;
	 // dwDataPos2 = dwDataPos1+dwDataHeight2;
	 // eax = &data[dwDataPos1].r
	 // ebx = &data[dwDataPos2].r
		 
	 mov   eax,dwDataHeight
	 mov   ebx,dwDataHeight2
	 imul  eax,edi
	 add   ebx,eax
	 
	 shl   eax,4  // eax*=16  COMPLEX 型は 64*2 bit変数
	 shl   ebx,4  // ebx*=16  
	 add	  eax,data  // ベースを足す
	 add	  ebx,data

L_HEIGHT:

	 // dReal1 = data[dwDataPos1].r;    // 実数 1
	 // dReal2 = data[dwDataPos2].r;    // 実数 2
	 // dImag1 = data[dwDataPos1].i;    // 虚数 1
	 // dImag2 = data[dwDataPos2].i;    // 虚数 2
	 
	fld qword ptr [eax]  // data[dwDataPos1].r
	fld qword ptr [ebx]  // data[dwDataPos2].r
	fld qword ptr [eax+8] // data[dwDataPos1].i
	fxch st(1)  // キャッシュを考えると fxch 入れた方が良い？
	fld qword ptr [ebx+8] // data[dwDataPos2].i

	// スタックの状態 st = i2,r2,i1,r1

	// data[dwDataPos1].i = dImag1+dImag2; // 虚数の和
	fld st(2)
	fadd st,st(1)
	fstp qword ptr [eax+8]  

	// data[dwDataPos1].r = dReal1+dReal2; // 実数の和
	fld st(3)
	fadd st,st(2)
	fstp qword ptr [eax]

	// st = i2,r2,i1,r1

	//虚数積
	// dReal3 = dReal1 - dReal2;
	// dImag3 = dImag1 - dImag2;

	fsubp st(2),st
	fsubp st(2),st

	// st = i1-i2,r1-r2 = img,real

	//  dCoss = dCosTable[dwStage*dwFftPoint/2+k];
	//  dSinn = dSinTable[dwStage*dwFftPoint/2+k];

	push ebx

	// ポインタ計算
	mov  ecx,esi
	mov  ebx,dwStage

	shl  ecx,3 // double 型なので *8 倍
	imul ebx,dwFftPoint

	shl  ebx,2

	add  ebx,ecx

	// Cos 値 Sin 値 FPU に代入
	mov  edx,dSinTable
	mov  ecx,dCosTable
	add edx,ebx
	add ecx,ebx

	pop ebx


	fld	qword ptr [edx]  //sin
	fld	qword ptr [ecx]  //cos
	fld st(1)
	fld st(1)

	// st = cos,sin,cos,sin,imag,real

  //  data[dwDataPos2].r = dReal3 * dCoss + dImag3 * dSinn;
  //  data[dwDataPos2].i = dImag3 * dCoss - dReal3 * dSinn;

	fmul st,st(5)
	fxch st(1)
	fmul st,st(4)
	faddp st(1),st
	fstp  qword ptr [ebx]

	// st = cos,sin,imag,real

	fmul st,st(2)
	fxch st(1)
	fmul st,st(3)
	fsubp st(1),st
	fstp  qword ptr [ebx+8]
	fstp st(0)
	fstp st(0)

	// st = 無し

	add eax,16	// eax = data[dwDataPos1] += 16
    add ebx,16	// ebx = data[dwDataPos2] += 16

	// ループ終わり
	inc esi
	cmp esi,dwDataHeight2
	jnae L_HEIGHT

	inc edi
	cmp edi,dwPart
	jnae L_PART

	// パートを 2 つに分ける
	// dwPart *=2; 
	// dwStage-- してループ
	mov ebx,dwPart
	mov eax,dwStage
	shl ebx,1
	dec eax
	mov dwPart,ebx
	mov dwStage,eax
	jnz L_STAGE
}

	//--------------------------------
	// 要素入れかえとパワースペクトル計算
	for(i = 0; i <  dwFftPoint ;i+=2){
		
		//ビットを反転して代入する位置を計算
		dwDataPos = 0;
		j = i;
		for(k=0;k< dwMaxBit;k++){
			dwDataPos << = 1;
			dwDataPos |= (j & 1);
			j >>= 1;
		}
		
		// パワースペクトルを計算して代入
		data[dwDataPos].r
			= data[i].r * data[i].r	+ data[i].i * data[i].i;
	}
	
}


//------------------------------------
// おまけ 
// パワースペクトル計算(C 版)
VOID CalcPowerSpec(COMPLEX* data,// FFT を行うデータ
				   DWORD dwFftPoint, // FFT の点数
				   DWORD dwMaxBit,  // dwFftPoint の最大ビット数
				   double* dCosTable,// Cos のテーブル
				   double* dSinTable// Sin のテーブル
				   ){

	// FFT 用変数
	DWORD dwStage;
	DWORD dwPart;
	DWORD dwDataHeight,dwDataHeight2;

	// データのポインタ
	DWORD dwDataPos1,dwDataPos2,dwPoint;

	// 実数、虚数
	double dReal1,dReal2,dReal3,dImag1,dImag2,dImag3;

	// 雑用
	DWORD i,j,k;
	double dSinn,dCoss;
		
	//------------
	// 計算部
	dwPart = 1;
	for(dwStage = (DWORD)(log(dwFftPoint)/log(2));dwStage>0;dwStage--)
	{
		dwDataHeight = (DWORD)1 <<  dwStage; // 2^dwStage
		dwDataHeight2 = dwDataHeight/2;
		
		//--------------------------------
		// バタフライ演算開始
		for(j = 0; j <  dwPart ;j++){	

			dwDataPos1 = dwDataHeight*j;
			dwDataPos2 = dwDataPos1+dwDataHeight2;

			for(k = 0 ; k <   dwDataHeight2 ; k++){
				
				dReal1 = data[dwDataPos1].r;	// 実数 1
				dImag1 = data[dwDataPos1].i;	// 虚数 1
				
				dReal2 = data[dwDataPos2].r;	// 実数 2
				dImag2 = data[dwDataPos2].i;	// 虚数 2
				
				// 足し合わせ	
				data[dwDataPos1].r = dReal1+dReal2;	// 実数の和
				data[dwDataPos1].i = dImag1+dImag2; // 虚数の和
				
				//虚数積			  
				dReal3 = dReal1 - dReal2;	  
				dImag3 = dImag1 - dImag2;  	   
				dwPoint = dwStage*dwFftPoint/2+k;
				dCoss = dCosTable[dwPoint];
				dSinn = dSinTable[dwPoint];
				
				data[dwDataPos2].r = dReal3 * dCoss + dImag3 * dSinn;	  
				data[dwDataPos2].i = dImag3 * dCoss - dReal3 * dSinn;

				dwDataPos1++;
				dwDataPos2++;
			}
		}
		
		dwPart *=2; // パートを 2 つに分ける
	}
	// 計算部ここまで
	//---------------

	//--------------------------------
	// 要素入れかえとパワースペクトル計算
	for(i = 0; i <  dwFftPoint ;i+=2){
		
		//ビットを反転して代入する位置を計算
		dwDataPos1 = 0;
		j = i;
		for(k=0;k< dwMaxBit;k++){
			dwDataPos1 << = 1;
			dwDataPos1 |= (j & 1);
			j >>= 1;
		}
		
		// パワースペクトルを計算して代入
		data[dwDataPos1].r
			= data[i].r * data[i].r	+ data[i].i * data[i].i;
	}
}