ビット数を数えるルーチン

32bitの立っているビット数を高速に数えるアルゴリズムとして、

int bitcount(unsigned long n)
{
    n = ((n&0xaaaaaaaa) >> 1) + (n&0x55555555); 
    n = ((n&0xcccccccc) >> 2) + (n&0x33333333); 
    n = ((n&0xf0f0f0f0) >> 4) + (n&0x0f0f0f0f); 
    n = ((n&0xff00ff00) >> 8) + (n&0x00ff00ff); 
    n = ((n&0xffff0000) >> 16) + (n&0x0000ffff); 
    return n;
}

こういうのが良く知られている。1行目で2ビットずつ足し合わせることで2ビットの中の立っているビット数を求め、2行目でそれらを2つずつ合計して4ビットの中の立っているビット数を求め…というのを32bit使って並列にやっている訳だ。ポイントを理解すれば単純明快、1ビットずつ32回評価するより段違いに速いのは明らかである。

数年前にこれを初めて見た時、独学の計算量節約プログラミングを長年続け、私の頭の中で完成しつつあった独自のプログラム理論は、一気に消し飛んだ。

そんな記憶も彼方の昨日、ネットサーフィンしていて、下のようなコードを見つけた。HAKMEM 169と呼ばれるアルゴリズムらしい。ネット上の色々なC言語実装を参考に勝手に整形した。

int bitcount2(unsigned long n){
    n = n - ((n >> 1) & 033333333333) - ((n >> 2) & 011111111111);
    return ((n + (n >> 3)) & 030707070707) % 077;
}

短い！
１つ目の例と同じ香りが漂うが、演算子の数が少ない。もしかして最強コードか？と思って解読しようとすると、数学的に大変ややこしい。
１行目は、3bitの整数m内の立っているビット数がm - (m>>1) - (m>>2)で計算できる（※１）ことを利用して、3bitずつ並列に、立っているビット数を計算する。２行目は、3bitずつのビット数を２つずつ足し合わせて6bitずつのビット数にしてから（※２）、bが6bitの整数、kが6の倍数の時にb*2^kを6bitの111111で割った余りがbになる（※３）ことを利用して、6bitずつの合計を求めている。

よく読むと、剰余演算を使ってるので、1つ目の例の方が速い。
参考：http://infolab.stanford.edu/~manku/bitcount/bitcount.html
（"Parallel"が１つ目の例、"MIT"が２つ目の例に相当する。"Parallel"より速いのはビット数のテーブルを予め用意している。）

元々のMITのHAKMEMのコードは剰余演算のあるアセンブラで書かれており、アセンブラとしては目を疑うほど短い。HAKMEM 169は最速ではないが、コードのコンパクトさと合わせて称えるべきアルゴリズムなのだろう。

※１：どのビットも、1つずつ単独で見ると、1つ右にシフトした値を引くと1つ右に移り、繰り返すと1の位に移り、それ以上は動かない。それを重ね合わせると、mが何ビットでも、m - (m>>1) - (m>>2) - (m>>3) ...が1のビット数になる。各ビットを1の位に移す数を同時に引くので、各ビットの移動が干渉せず、各ビットを1の位に移したものの和になる。

※２：1のビット数の最大が32なので、結果は6bitの数であり、3bitに区切ったままでは一気に合計を出せないため、6bit区切りにする。

※３： yを自然数、xを0<=x<2^yの整数とすると、x*2^y - x*(2^y-1) = x, x*2^{2y} - x*2^y*(2^y-1) = x*2^yであることから、x*2^{yの倍数}を(2^y-1)で割った余りがxになることがわかる。

※なお、ここに英語ながら割とわかりやすい説明があるが、%63する対象の説明は間違っており、読むと混乱する。（同ページのコメント欄参照）