基本的な演算(加算、減算、論理演算、シフト) を組み合わせて 応用的な演算(乗除算・ログ・二乗、ビット操作、バイト操作、浮動小数演算) を いかに実現するかというアルゴリズムの本。
でも、どのアルゴリズムも超絶テクニックばっかりです。
前文はハッカー界の界王様 Guy L. Steele, Jr が寄せている。
作成日：2004.05.04
修正日：2006.03.18
このページは 2003年の 9/11、 9/28 の日記をまとめて作成。
PowerPC 系や Alpha などには population count と呼ばれる
レジスタ中の立っているビット数を数える命令が実装されている。
集合演算を行うライブラリを実装したい場合などに重宝しそうな命令である。
職場でこの population count 命令について話をしているうちに
ビットカウント操作をハードウェアで実装するのは得なのか？という点が議論になった。
CPU の設計をできるだけシンプルにするためには、
複雑で使用頻度の低い命令は極力減らした方がよい。
例えば SPARC は命令セット中にビットカウント演算があるが、
CPU 内には実装しないという方針をとっている
(population 命令を実行すると不正命令例外が発生し、
それを OS がトラップして処理する)。
だが管理人はビットカウント演算が十分に使用頻度が高い場合には
CPU がハードウェア処理した方がペイするのではないか？と思っていた。
しかし識者の意見では、 ビットカウント演算はビット演算命令の組み合わせで ハードウェア回路を用意した場合 以上の速度が出せるのでまったく無駄であるとのこと。 実際にビット数を数えるアルゴリズムのコードを見てかなりショックを受けた。
このページにそのアルゴリズムを記す。 特に断らない場合、プログラムは 32-bit CPU での実行を念頭においている。
レジスタ中の 1 になっているビット数を数えるアルゴリズムについて、
バージョン 1 から 5 までを紹介する。
後になるほど最適化されて行くのが分かると思う。
整数変数中のビットをカウントするプログラムとしては、
プログラマーなら誰でも
以下のようなプログラムが反射的に考えつくはず。
ただ
このアルゴリズムは非常に大きな時間コストが掛かる
ダメアルゴリズムだ。
int numofbits1(int bits) {
int num = 0;
int mask = 1;
for( ; mask != 0 ; mask = mask << 1 ){
if( bits & mask )
num++ ;
}
return num;
}
次に考えつくのは、
バージョン2のようなテーブルを用いたやり方だろう。
実際にはあまり大きなテーブルは作成できないので
0 〜 255 までの 8 ビット分のビットカウントテーブルを作成しておいて、
ビットカウントしたデータを 8 ビットづつに区切って
処理することになるだろう。
(ループをまわす回数が一定ならループを展開しても構わないので、
条件分岐は不要となるかもしれない。)
const int BITS_COUNT_TABLE[256] = { /* 予めプリセット */ } ;
int numofbits2(int bits) {
int num = 0;
for( int i=0 ; i<sizeof(bits) ; i++ ){
num += BITS_COUNT_TABLE[((unsigned char*)&bits)[i]]
}
return num;
}
以下のようなアルゴリズムまでは
自力で考えつくかも知れない。
このアルゴリズムでは
立っているビットが疎な場合には
かなり素早く収束する。
int numofbits3(int bits) {
int num = 0;
for( ; bits != 0 ; bits &= bits - 1 ){
num++ ;
}
return num;
}
ここからが管理人がショックを受けたアルゴリズム。
バージョン 4 ではループを使用しないでも
ビットカウント操作ができてしまっている。
int numofbits4(int bits) {
int num;
num = (bits >> 1) & 03333333333;
num = bits - num - ((num >> 1) & 03333333333);
num = ((num + (num >> 3)) & 0707070707) % 077;
return num;
}
ビットを数える最適化されたアルゴリズム。
このアルゴリズムではループ・条件分岐がないばかりか、
剰余命令まで消失している。
現在のスーパースカラープロセッサでは
演算を畳み込めるので、
CPU にビルトインしたビットカウント命令よりも
このアルゴリズムを用いた方が
高速に処理ができると思われる。
int numofbits5(long bits) {
bits = (bits & 0x55555555) + (bits >> 1 & 0x55555555);
bits = (bits & 0x33333333) + (bits >> 2 & 0x33333333);
bits = (bits & 0x0f0f0f0f) + (bits >> 4 & 0x0f0f0f0f);
bits = (bits & 0x00ff00ff) + (bits >> 8 & 0x00ff00ff);
return (bits & 0x0000ffff) + (bits >>16 & 0x0000ffff);
}
上記のアルゴリズムは、 バージョン4 が Knuth の本に、 バージョン5 が 参考文献 にあげた Hacker's Delight に載っている。
人類がバージョン5 までたどり着いたのは 1960 年代だったようだ。
IBM の研究者が
IBM Stretch 用のプロセッサのために
count population や count leading zero に関して
徹底的な調査を試みた古典的な論文があるそうだ
(Count leading zero は、
ビット列を MSB から順に読んで最初に1になったビットが現れるまで
0 がいくつ連続して出てくるかを数える操作。
LSB から数えるのは count trailing zero)。
当然 count leading zero にも高速アルゴリズムがある。
悔しいので答を見ずに考えてみることにする。
レジスタの中のビット列の左側 Leftmost bit から見て 0 が連続している数を Number of Leading Zero (NLZ) という (0 の場合には 32 が返る)。 言い換えると 0 ビット目から見て 1 が最初に立っているビット位置のことである (1 がどこにも立っていない場合には 32 が返るとする)。 同様にレジスタの中のビット列の右側 Rightmost bit からみて 0 がいくつ連続しているかを Number of Training Zero (NTZ) と呼ぶ。
例をあげると
|
32-bit 値 0000010010000100 の場合、 NLZ は 5 で、NTZ は 2 となる。
NTZ は比較的簡単である。
x&(-x) を行うと左側から最初に 1 が立っているアドレスのみを残して
他のビットは 0 で埋まるので、
(x&(-x))-1 を
計算しそのビット数を数えればよい。
|
NTZ は NTZ(x) = count_bits((~x)&(x-1)) でもよい。
NLZ に関してはあまりスマートなアルゴリズムはないが、 3 つ (+1) のアルゴリズムを挙げる。
まず、 すでに最適化な CLZ がすでに見つかっている場合には、
NTZ(x) = 32 - NLZ((~x)&(x-1))
がなりたつ。 そのため、 以下のような方法で
int nlz1(unsigned x) {
x = x | ( x >> 1 );
x = x | ( x >> 2 );
x = x | ( x >> 4 );
x = x | ( x >> 8 );
x = x | ( x >> 16 );
return count_bits( ~x );
}
nlz1 の ~x までの操作で、
上位のビットが下位にコピーされて行くので
以下のようにビットが立つことになる。
|
入力 出力 0000 → 1111 0001 → 1110 0010 → 1100 0011 → 1100 0100 → 1000 0101 → 1000 0110 → 1000 0111 → 1000 1000 → 0000 1001 → 0000 |
ただし nlz1 は ほとんどのアーキテクチャで最速ではないようだ。
int nlz2(unsigned x) {
unsigned y;
int n = 32;
y = x >> 16; if (y != 0){ n = n - 16 ; x = y; }
y = x >> 8; if (y != 0){ n = n - 8 ; x = y; }
y = x >> 4; if (y != 0){ n = n - 4 ; x = y; }
y = x >> 2; if (y != 0){ n = n - 2 ; x = y; }
y = x >> 1; if (y != 0){ return n - 2; }
return n - x;
}
int nlz3(unsigned x) {
int y, m, n;
y = - (x >> 16);
m = (y >> 16) & 16;
n = 16 - m;
x = x >> m;
y = x - 0x100;
m = (y >> 16) & 8;
n = n + m;
x = x << m;
y = x - 0x1000;
m = (y >> 16) & 4;
n = n + m;
x = x << m;
y = x - 0x4000;
m = (y >> 16) & 2;
n = n + m;
x = x << m;
y = x >> 14;
m = y & ~(y >> 1);
return n + 2 - m;
}
NLZ を求める 3 つのアルゴリズムを Linux の AMD Duron 800MHz / gcc 2.95.3 (-O2) 環境で 実行し実測してみた。 すると
nlz2(7.7秒) >nlz3(11.7秒) >nlz1(18.1秒)
と分岐命令の多い nlz2 が最速という結果になった
(括弧内は 0x10000000 回 実行した実行時間)。
無論、結果はコンパイラやアーキテクチャによって
左右されると思われる
(分岐ペナルティが大きい Pentium4 では nlz2 は不利？)。
NLZ を求めるには 浮動小数点演算を利用する方法もある。
int nlz4(unsigned int x) {
int n;
#if USE_DOUBLE
union {
unsigned long long as_uint64;
double as_double;
} data;
data.as_double = (double)x + 0.5;
n = 1054 - (data.as_uint64 >> 52);
#else
union {
unsigned int as_uint32;
float as_float;
} data;
data.as_float = (float)x + 0.5;
n = 158 - (data.as_uint32 >> 23);
#endif
return n;
}
これは IEEE 754 形式の浮動小数点フォーマットを利用した方法。
倍精度浮動小数(64ビット)のメモリイメージは
以下のようなビット配置になり、
その値は (-1)^S × 1.xxxxxx × 2^(yyyy) の形になるので、
整数を浮動小数点に変換すると一番左に 1 が立っている位置が仮数部に現れる。
| 63 | 62〜52 | 51〜0 |
| 符号(S) | 指数(yyy...) | 仮数(xxx...) |
これは結局、整数 → 浮動小数への変換が
log2 の演算に相当するから。
nlz4 で 0.5 を足すのは x = 0.0 の時
正しく 32 を返させるためのおまじない。
浮動小数点の 0.0 のメモリ表現は
すべてのビットが 0 となる特殊な値なので、
1.0 × 2^(-1) となる
0.5 をストッパーとして足しておく。
整数 浮動小数点 0001 → 1.000 × 2^(0) 001x → 1.x00 × 2^(1) 01xx → 1.xx0 × 2^(2) 1xxx → 1.xxx × 2^(3)
ただし C 版の nlz4 のコードを実際にコンパイルすると、
共用体部分に load / store 命令部分が出現するため
あまり速度は出てないと予想される。
SPARC 系を除く大概の CPU には
整数レジスタ ⇔ 浮動小数点レジスタ命令があるので、
アセンブラを使ってごりごり記述するが吉。
あとこの方法はかなり奇天烈だが、 結局のところビットの検索を 浮動小数点演算器というブラックボックスに任せただけなので アルゴリズムとは言い難いかも。
基本的な演算(加算、減算、論理演算、シフト) を組み合わせて 応用的な演算(乗除算・ログ・二乗、ビット操作、バイト操作、浮動小数演算) を いかに実現するかというアルゴリズムの本。
でも、どのアルゴリズムも超絶テクニックばっかりです。
前文はハッカー界の界王様 Guy L. Steele, Jr が寄せている。
http://www.amazon.co.jp/exec/obidos/ASIN/4434046683/250-9184073-6450663
nminoruさんには悪いけど原著をもう買ったので、多分買わないだろうなぁ。
会社の資料室には買わせるかも。
邦訳は 3,570円 と価格的にお徳なので、原著を持っていない人にはよいかも (^_^;