多倍長整数でSteinのGCDを使う
| プログラムメモ | 2012/04/22 Sun 22:43
次のような、32bit系の多倍長整数クラスがあったとして、
16進で次のように表せる整数があれば
0x111100002222000033330000
このクラスでは次のように格納されるものとする。(注1)
これを用いて、Stein's algorithmによるGCD計算を考える。
SteinのGCDについてはWikiPediaが詳しい。
http://en.wikipedia.org/wiki/Binary_GCD_algorithm
一般に ユークリッドの互除法(拡張) の方がループ回数は少ないが、
剰余計算がネックになり、応答性が悪くなる場合が多い。
但し、実装する乗算・剰余算のアルゴリズムや数値の桁数によっては、
ループ回数の少ないユークリッドの互除法(拡張)の方が早い場合もある。
Steinのアルゴリズムで遅いと思われるのは、Shift演算の回数に尽きる。
サンプルコードをそのまま使うと、いくら負荷の少ないShift演算と言えども、
多桁演算で処理回数が多くなると大きな遅延を招きます。
WikiPediaのSteinのサンプルコードで問題となるのは、
の部分。要は右端の０を除去しているだけである。
右端に並んでいる０の数を数える関数を用意して、Shift回数を減らことを目指す。
http://en.wikipedia.org/wiki/Count_trailing_zeros にもあるが、
intel系であれば、アセンブラ(bsf)が使える。
BITPERBLOCKは、配列１つに含まれる多倍長整数のBit数である。
32bitであれば、32が使われる場合が多い。(注2)
(64bitの場合、bsfqを使う。)
(m_sz=2, BITPERBLOCK=32)の次のような多倍長整数があるとき
アセンブラを使わない次の方法もある。
http://en.wikipedia.org/wiki/Count_trailing_zeros にもありますが、
Ｍ系列を使ったマジックです。
こちらの方が詳しいです。 http://d.hatena.ne.jp/siokoshou/20090704
すると、Steinのコードは次のように書ける。
これで30〜40%は、処理時間を短縮できました。
これこそがSteins;Gateの選択♪
このコードを実装したサンプルの16進・２進・10進電卓です。
整数だけですがおよそ100万桁まで計算できます。
(注1) 固定長ではないBigintクラスでは、伸長・縮退等のメモリの管理も必要になるので、これほど単純なものではありません。
(注2) 32bit(Block)全てを使うことが必ずしも有利ではありません。SSEを有効に使うなら桁あふれが発生しないよう 31bit 以下にします。
Tags: プログラムメモ
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)の次のような多倍長整数があるとき
bsf()は、32+16=48 を返す。|- m_num[0]-|- m_num[1]-| | 0000 0000 | ABCD 0000 |
アセンブラを使わない次の方法もある。
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: プログラムメモ
author : HUNDREDSOFT | - | -
