階乗の計算入門から精通へ―入門編の三組の威力

要約:本稿では、乗算を計算するプログラムを2つ提供します.1番目のプログラムはCにアセンブラコードを埋め込む方法を採用し,前編でプログラム2のボトルネック部分を改善し,速度を従来の3倍以上に向上させた.2番目のプログラムはアルゴリズムをさらに改良し、1万の階乗を計算する際に、前編のプログラムより2-6倍速く、10000の階乗を計算し、
1.7 Gを疾走するには0.25秒しかかかりません.
     
前の文章で「大数階乗の計算は入門から精通まで」
-入門編の2』では,64ビット整数の除算を用いるため,2番目のプログラムはかえって速度が遅い2つの大きな数次乗を計算するプログラムを与えた.本稿では,C言語にアセンブリコードを埋め込む手法を用いて,ボトルネック部分を改善し,計算速度を従来の3倍以上に向上させた.
まず,大数次乗を計算するコアコード(以下を参照)を見ると,ループごとに1回64ビットの乗算,1回64ビットの加算,2回64ビットの除算が必要であることがわかる.
余を求めれば除算と見なす).除法命令は乗算命令よりも遅く,64ビットの除法ではなおさらであることが分かった.にある
vcでは、2つの64ビット数の除算は変調です.
aulldiv関数で実現された2つの64ビット数
の余剰演算は呼び出しです
aullremで実現しました.この2つの関数のソースコード(アセンブリコード)を解析することにより、除算演算を行う際に、まず除算数をチェックし、2の32次方未満であれば、64ビット商を得るために2回の除算が必要であり、除算数が2以上であれば
32、演算はさらに複雑になります.同じように
余剰演算をする場合、除数が2未満の場合
32,32ビットの残数を1つ得るためにも2回の除算が必要である.この例では、除算数は10000000で、2未満です.
32です.そのため、このコードは4回の除算が必要です.
 

for (carry=0,p=pTail;p>=pHead;p--)
{
    prod=(UINT64)(*p) * (UINT64)i +carry;
    *p=(DWORD)(prod % TEN9);
    carry=prod / TEN9;
}

このコードでは,除算演算を行う場合,その商は常に2の32乗未満であるため,この64ビット数の除算は
および余剰演算は、1つの除算命令で実現することができる.次に,C関数にアセンブリコードを埋め込む方法で,1回の除算命令を同時に実行する.
商数を得る.関数の原形は次のとおりです.
DWORD Div_TEN9_2 (UINT64 x,DWORD *pRemainder );これ
関数はxを10000000で割った商を返して、除数は預け入れます
pRemainder.次は関数Div_TEN9_2と階乗を計算するコードは,Cに埋め込まれたアセンブリで実現される.Cコードに埋め込まれたアセンブリの使い方についてはMSDNを参照してください.
  

inline DWORD Div_TEN9_2(UINT64 x,DWORD *pRemainder )
{
	_asm
	{
		mov eax,dword ptr [x]   // low DWORD  
			mov edx,dword ptr [x+4] // high DWORD  
			
			mov ebx,TEN9
			div ebx
			mov ebx,pRemainder
			mov [ebx],edx                 // remainder-> *remainder 
			// eax, return value 
	}
}

void calcFac2(DWORD n)
{
	DWORD *buff,*pHead,*pTail,*p;
	DWORD t,i,len,carry;
	UINT64 prod;
	
	if (n==0)
	{ printf("%d!=1",n); return; }
	
	//---             
	t=GetTickCount();
	
	len=calcResultLen(n,TEN9);
	buff=(DWORD*)malloc( sizeof(DWORD)*len);
	if (buff==NULL)
		return ;
	
	//      n!  
	pHead=pTail=buff+len-1;
	for (*pTail=1,i=2;i<=n;i++)
	{
		for (carry=0,p=pTail;p>=pHead;p--)
		{
			prod=(UINT64)(*p) * (UINT64)i +(UINT64)carry;
			carry=Div_TEN9_2(prod,p );
		}
		
		while (carry>0)
		{
			pHead--;
			*pHead=(DWORD)(carry % TEN9);
			carry /= TEN9;
		}
	}
	
	t=GetTickCount()-t;
	
	//       
	printf("It take %d ms/n",t);
	printf("%d!=%d",n,*pHead);
	for (p=pHead+1;p<=pTail;p++)
		printf("%09d",*p);
	printf("/n");
	
	free(buff);    //        
}

なお,本稿の題名はアセンブリの威力であるが,アセンブリ言語では必ずしも有効ではないが,一般的にはキーコードをアセンブリで書き換えた場合,性能向上の幅は上記の例ほど顕著ではなく,せいぜい30%向上する可能性があり,本プログラムは特例である.
 
最後の改良:これらの計算階乗の関数を解析すると、計算階乗は実際には二重サイクルであり、内サイクル部分は(i-1)!*
i,外ループは順次2を計算する!
,3!,4!,nまで!もし私たちがr=(i-1)!,prod=を先に算出できるかどうか
i*(i+1)*...mは、
i*(i+1)*...ちょうど2^32未満で、
i*(i+1)*...m*(m+1)では>=2^32となり、r*prodを再計算することで、外部サイクルの回数を減らし、速度を向上させることができる.理論的および試験的結果は,30000以下の乗算を計算すると,速度が1倍以上向上できることを示した.コードを次に示します.
 

void calcFac3(DWORD n)
{
	DWORD *buff,*pHead,*pTail,*p;
	DWORD t,i,len,carry;
	UINT64 prod;
	
	if (n==0)
	{ printf("%d!=1",n); return; }
	
	//---             
	t=GetTickCount();
	
	len=calcResultLen(n,TEN9);
	buff=(DWORD*)malloc( sizeof(DWORD)*len);
	if (buff==NULL)
		return ;
	
	//      n!  
	pHead=pTail=buff+len-1;
	for (*pTail=1,i=2;i+15<n;)
	{
		UINT64 t=i++;
		while (t<4294967296I64)
		{
			t *= (UINT64)i;                 i++;
		}
		i--;
		t/=i;
		
		for (carry=0,p=pTail;p>=pHead;p--)
		{
			prod=(UINT64)(*p) * t +(UINT64)carry;
			carry=Div_TEN9_2(prod,p );
		}
		
		while (carry>0)
		{
			pHead--;
			*pHead=(DWORD)(carry % TEN9);
			carry /= TEN9;
		}
	}
	
	for (;i<=n;i++)
	{
		for (carry=0,p=pTail;p>=pHead;p--)
		{
			prod=(UINT64)(*p) * (UINT64)i +(UINT64)carry;
			carry=Div_TEN9_2(prod,p );
		}
		
		while (carry>0)
		{
			pHead--;
			*pHead=(DWORD)(carry % TEN9);
			carry /= TEN9;
		}
	}
	
	printf("It take %d ms/n", GetTickCount()-t);
	//       
	printf("%d!=%d",n,*pHead);
	for (p=pHead+1;p<=pTail;p++)
		printf("%09d",*p);
	printf("/n");
	
	free(buff);//        
}

[email protected]、著作権所有、転載は出典を明記してください.