Hundredsoft

次のような、32bit系の多倍長整数クラスがあったとして、

class Bigint
{
public:
	Bigint();
	...
	... 様々な定義(四則演算、operator=,operator<<=等)
	...
private:
	int   m_sz;   // 配列数
	int*  m_num;  // 配列(リトルエンディアン)
};

16進で次のように表せる整数があれば
0x111100002222000033330000

このクラスでは次のように格納されるものとする。(注1)

(m_sz=3)
|- m_num[0]-|- m_num[1]-|- m_num[2]-|
| 3333 0000 | 2222 0000 | 1111 0000 |

これを用いて、Stein's algorithmによるGCD計算を考える。


SteinのGCDについてはWikiPediaが詳しい。
http://en.wikipedia.org/wiki/Binary_GCD_algorithm

一般に ユークリッドの互除法(拡張) の方がループ回数は少ないが、
剰余計算がネックになり、応答性が悪くなる場合が多い。

但し、実装する乗算・剰余算のアルゴリズムや数値の桁数によっては、
ループ回数の少ないユークリッドの互除法(拡張)の方が早い場合もある。

Steinのアルゴリズムで遅いと思われるのは、Shift演算の回数に尽きる。
サンプルコードをそのまま使うと、いくら負荷の少ないShift演算と言えども、
多桁演算で処理回数が多くなると大きな遅延を招きます。
WikiPediaのSteinのサンプルコードで問題となるのは、

/* From here on, u is always odd. */
do {
  while ((v & 1) == 0)  /* Loop X */
     v >>= 1;

の部分。要は右端の０を除去しているだけである。
右端に並んでいる０の数を数える関数を用意して、Shift回数を減らことを目指す。

http://en.wikipedia.org/wiki/Count_trailing_zeros にもあるが、
intel系であれば、アセンブラ(bsf)が使える。

//
// 高速化のためゼロチェックを行っていない
//  ゼロだとアクセス違反の例外。
//
int Bigint::bsf() const{
	int* p = m_num;
	__asm{
		mov edx, [p]
		mov eax, edx
bsf_l2:
		bsf ecx, [eax]
		jnz bsf_l1
		add eax, 4
		jmp bsf_l2
bsf_l1:
		sub eax, edx
		shr eax, 2
		mov edx, BITPERBLOCK
		mul edx
		add eax, ecx
	}
}

BITPERBLOCKは、配列１つに含まれる多倍長整数のBit数である。
32bitであれば、32が使われる場合が多い。(注2)
(64bitの場合、bsfqを使う。)

(m_sz=2, BITPERBLOCK=32)の次のような多倍長整数があるとき

|- m_num[0]-|- m_num[1]-|
| 0000 0000 | ABCD 0000 |

bsf()は、32+16=48 を返す。


アセンブラを使わない次の方法もある。

int Bigint::bsf() const{
	static int table[32] = {0,1,10,2,11,14,22,3,30,12,15,17,19,23,26,4,31,9,13,21,29,16,18,25,8,20,28,24,7,27,6,5};

	int count = 0;
	for (int i=0; i<m_sz; i++){
		int x = m_num[i];
		if (x == 0){
			count += BITPERBLOCK;
		}else{
			count += table[((unsigned int)(x & -x) * 0x07C4ACDDUL) >> 27];
			break;
		}
	}
	return count;
}

// 64bit/Blockの場合
int Bigint::bsf64() const{
	static int table[64] = {
		0,1,59,2,60,40,54,3,61,32,49,41,55,19,35,4,62,52,30,33,50,12,14,42,56,16,27,20,36,23,44,5,
		63,58,39,53,31,48,18,34,51,29,11,13,15,26,22,43,57,38,47,17,28,10,25,21,37,46,9,24,45,8,7,6};

	int count = 0;
	for (int i=0; i<m_sz; i++){
		__int64 x = m_num[i];
		if (x == 0){
			count += BITPERBLOCK;
		}else{
			count += table[((unsigned __int64)(x & -x) * 0x03F566ED27179461ULL) >> 58];
			break;
		}
	}
	return count;
}

http://en.wikipedia.org/wiki/Count_trailing_zeros にもありますが、
Ｍ系列を使ったマジックです。
こちらの方が詳しいです。 http://d.hatena.ne.jp/siokoshou/20090704


すると、Steinのコードは次のように書ける。

//
// *this = GCD(*this, v)
//
void Bigint::gcd(const Bigint& v) {
	if (v == 0){
		return;
	}
	if (operator==(0)){
		operator= (v);
		return;
	}

	Bigint x[2] = {*this, v};

	if (x[0] < 0)	x[0].operator-();
	if (x[1] < 0)	x[1].operator-();

	int i0 = x[0].bsf();
	int i1 = x[1].bsf();
	int shift = (i0 < i1) ? i0 : i1;
	x[0] >>= i0;
	x[1] >>= shift;

	int id = 0;
	do {
		x[id^1] >>= x[id^1].bsf();

		if (x[id] > x[id^1]){
			id ^= 1;
		}
		x[id^1] -= x[id];
	} while (x[id^1] != 0);

	operator= (x[id] << shift);
}

これで30〜40%は、処理時間を短縮できました。

これこそがSteins;Gateの選択♪

このコードを実装したサンプルの16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。


(注1) 固定長ではないBigintクラスでは、伸長・縮退等のメモリの管理も必要になるので、これほど単純なものではありません。
(注2) 32bit(Block)全てを使うことが必ずしも有利ではありません。SSEを有効に使うなら桁あふれが発生しないよう 31bit 以下にします。

---

Tags: プログラムメモ