ビットを数える・探すアルゴリズム

作成日：2004.05.04
修正日：2012.09.01

このページは 2003年の9/11、9/28 の日記をまとめて作成。

はじめに

PowerPC 系や Alpha などには population count と呼ばれるレジスタ中の立っているビット数を数える命令が実装されている。 集合演算を行うライブラリを実装したい場合などに重宝しそうな命令である。

職場でこの population count 命令について話をしているうちにビットカウント操作をハードウェアで実装するのは得なのか？という点が議論になった。 CPU の設計をできるだけシンプルにするためには、複雑で使用頻度の低い命令は極力減らした方がよい。 例えば SPARC は命令セット中にビットカウント演算があるが、CPU 内には実装しないという方針をとっている(population 命令を実行すると不正命令例外が発生し、それを OS がトラップして処理する)。 だが管理人はビットカウント演算が十分に使用頻度が高い場合には CPU がハードウェア処理した方がペイするのではないか？と思っていた。

しかし識者の意見では、ビットカウント演算はビット演算命令の組み合わせでハードウェア回路を用意した場合以上の速度が出せるのでまったく無駄であるとのこと。 実際にビット数を数えるアルゴリズムのコードを見てかなりショックを受けた。

このページにそのアルゴリズムを記す。 特に断らない場合、プログラムは 32-bit CPU での実行を念頭においている。

1. ビット数を数えるアルゴリズム (Counting 1-bits)
2. MSB から (LSB から) 0 が立っている数を数える
3. 参考文献

1. ビット数を数えるアルゴリズム (Counting 1-bits)

レジスタ中の 1 になっているビット数を数えるアルゴリズムについて、バージョン 1 から 5 までを紹介する。
後になるほど最適化されて行くのが分かると思う。

バージョン1

整数変数中のビットをカウントするプログラムとしては、プログラマーなら誰でも以下のようなプログラムが反射的に考えつくはず。
ただこのアルゴリズムは非常に大きな時間コストが掛かるダメアルゴリズムだ。

int numofbits1(int bits)
{
    int num  = 0;
    int mask = 1;

    for ( ; mask != 0 ; mask = mask << 1 ){
        if (bits & mask )
            num++;
    }

    return num;
}

バージョン2

次に考えつくのは、バージョン2のようなテーブルを用いたやり方だろう。
実際にはあまり大きなテーブルは作成できないので 0 〜 255 までの 8 ビット分のビットカウントテーブルを作成しておいて、ビットカウントしたデータを 8 ビットづつに区切って処理することになるだろう。
(ループをまわす回数が一定ならループを展開しても構わないので、条件分岐は不要となるかもしれない。)

const int BITS_COUNT_TABLE[256] = {
    0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
    4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8,
};

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;
}

バージョン3

以下のようなアルゴリズムまでは自力で考えつくかも知れない。
このアルゴリズムでは立っているビットが疎な場合にはかなり素早く収束する。

int numofbits3(int bits)
{
    int num  = 0;

    for ( ; bits != 0 ; bits &= bits - 1 ) {
        num++;
    }

    return num;
}

バージョン4

ここからが管理人がショックを受けたアルゴリズム。
バージョン 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;
}

バージョン5

ビットを数える最適化されたアルゴリズム。
このアルゴリズムではループ・条件分岐がないばかりか、剰余命令まで消失している。 現在のスーパースカラープロセッサでは演算を畳み込めるので、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);
}

性能測定

バージョン1から5までと SSE 4.2 で追加された POPCNT 命令の性能を計測してみる。

• numofbits_routines.c
• numofbits_test.c

これを gcc 4.4.6 を使い以下のようにコンパイルする。

$ gcc -g -Wall -mtune=core2 -O3 numofbits_test.c numofbits_routines.c

Core i5-2400 3.1GHz マシンで性能を測定したところ、以下の結果が得られた。

numofbits1	27.376090 sec
numofbits2	0.657061 sec
numofbits3	16.755048 sec
numofbits4	2.434369 sec
numofbits5	1.678582 sec
popcnt	0.338981 sec

バージョン1 → 3 → 4 → 5 と性能がよくなっているのが分かる。 ただしテーブルを引くバージョン 2 は意外に性能がよい。

Intel x86 アーキテクチャが SSE 4.2 から導入した population count 命令の POPCNT は最速である。 良好なアルゴリズムがあってもハードウェアで実際してしまった方がだいぶ性能がよいようだ。

上記のアルゴリズムは、バージョン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 にも高速アルゴリズムがある。 悔しいので答を見ずに考えてみることにする。

MSB から (LSB から) 0 が立っている数を数える

レジスタの中のビット列の左側 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 を求める

NTZ は比較的簡単である。 x & (-x) を行うと左側右側から最初に 1 が立っているアドレスのみを残して他のビットは 0 で埋まるので、(x & (-x))-1 を計算しそのビット数を数えればよい。

NTZ は NTZ(x) = count_bits((~x) & (x-1)) でもよい。

NLZ を求める

NLZ に関してはあまりスマートなアルゴリズムはないが、3つ (+1) のアルゴリズムを挙げる。

バージョン1: ビットカウント問題に帰着させる場合

まず、すでに最適化な CLZ がすでに見つかっている場合には、

NTZ(x) = 32 - NLZ((~x) & (x-1))

がなりたつ。 そのため、以下のような方法で

int nlz1(unsigned int 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 はほとんどのアーキテクチャで最速ではないようだ。

バージョン2： 分岐を使う

int nlz2(unsigned int x)
{
    unsigned int 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;
}

バージョン3： 分岐を使わない

int nlz3(unsigned int 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;
}

バージョン4： 浮動小数点演算を利用

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 には整数レジスタ ⇔ 浮動小数点レジスタ命令があるので、アセンブラを使ってごりごり記述するが吉。

あとこの方法はかなり奇天烈だが、結局のところビットの検索を浮動小数点演算器というブラックボックスに任せただけなのでアルゴリズムとは言い難いかも。

性能測定

いろいろなバージョンの NLZ の性能を計測してみる。 測定プログラムは以下の通りで、

• nlz_routines.c
• nlz_test.c

これを gcc 4.4.6 を使い以下のようにコンパイルする。

$ gcc -g -Wall -O3 -mtune=core2 nlz_test.c nlz_routines.c

Core i5-2400 3.1GHz マシンで性能を測定したところ、以下の結果が得られた。

nlz1	ビットカウント問題に帰着	1.032112 sec
nlz2	分岐を使う	8.323925 sec
nlz3	分岐を使わない	4.341245 sec
nlz4	浮動小数点演算	0.352829 sec
nlz5	ささ氏の方法	1.587580 sec

最速はやはり浮動小数点演算器を使った nlz4 である。
nlz1 は SSE 4.2 の POPCNT 命令を使っているため非常に高速になっている。

特殊なハードを使わないアルゴリズムでは、コメント[7]のささ氏の手法が一番高速な結果が得られている。

参考文献

Henry S. Warren Hacker's Delight (2nd Edt. (1st Edt. の日本語訳 ハッカーのたのしみ)

基本的な演算(加算、減算、論理演算、シフト) を組み合わせて応用的な演算(乗除算・ログ・二乗、ビット操作、バイト操作、浮動小数演算) をいかに実現するかというアルゴリズムの本。
でも、どのアルゴリズムも超絶テクニックばっかりです。

前文はハッカー界の界王様 Guy L. Steele, Jr が寄せている。

コメント