一番右端の立っているビット位置を求める「ものすごい」コード

一番右端の立っているビット位置(RightMostBit)を求めるコードで速いのないかなーと探していたら、ものっっっすごいコードに出会ってしまったのでご紹介。2ch のビット演算スレで 32bit 値のコードに出会って衝撃を受けて、その後 64bit 値版のヒントを見つけたのでコードを書いてみました。
この問題は ハッカーのたのしみ―本物のプログラマはいかにして問題を解くか (Google book search で原著 Hacker's delight が読めたのでそれで済ませた) で number of trailing zeros (ntz) として紹介されています。bit で考えたときに右側に 0 がいくつあるかを数えるもの。1 だと 0、2 だと 1、0x80 なら 7、12 なら 2 といったぐあい。0 のときに表題どおりの問題として考えるといくつを返すの?ってことになるので、問題を正確にするために ntz として 64 を返すことにしたんだと思います。すぐに思いつくのがループで 1bit ずつ見ていく方法ですよね。高速化する方法はおもいつきますか?Hacker's delight ではバイナリーサーチを使ってて、見たときはスゲーと思ったんだけど、このコードの前では色あせて見えます(^^;
ちなみに環境依存してよいなら x86/x64 に bsf (VC++ なら http://msdn.microsoft.com/en-us/library/wfd9z0bb.aspx) ってのがあります。bsf では 0 のときの戻り値は未定義です。

問題の説明はここまでにして、コードの紹介です。Hacker's delight のコードより4〜5倍速く(5-13より4〜5倍、5-15より1.2〜1.3倍)、そして、イミフ加減が半端じゃない!これ一つで 64bit 値以下のすべての値に対応できます。

public static int GetNumberOfTrailingZeros( long x )
{
    if ( x == 0 ) return 64;

    ulong y = ( ulong ) ( x & -x );
    int i = ( int ) ( ( y * 0x03F566ED27179461UL ) >> 58 );
    return table[ i ];
}

table は以下のようにして求めます。

static int[] table;

table = new int[ 64 ];
ulong hash = 0x03F566ED27179461UL;
for ( int i = 0; i < 64; i++ )
{
    table[ hash >> 58 ] = i;
    hash <<= 1;
}

意味わかりますか?w
x & -x は Hacker's delight にある、立っている一番右端のビットだけ残して0にしてしまう黒魔術。使える場面が多いので覚えておくと便利です。いろんな値を入れてためしてみてください。
その後、完全ハッシュを使って数(2^n)を数(0〜63)に変換してるわけです。0のときの if 文が目障りですが、0のとき呼ばないなど don't care であれば省いてもOKです。省くと 0 のときに 0 が返ってきます。bit1 最下位ビットが立っているときも 0 なので区別できません。64bit 値の右端のビット位置をあらわすには、0の場合も含めると65通りの答えが必要で、もし65をあまり越えないコンパクトなハッシュ値を作れると if が不要になります。あってもこっちのほうが速い可能性はありますが。

このコードは 2ch のビット演算スレ0x03 の 71 と 80 で知りました。
http://pc12.2ch.net/test/read.cgi/tech/1226143920/

80 から解説の引用の引用。

周期2^p-1ビットの列があって、そこからpビットを切り出すと、オール0以外のすべてのパターンが現れる

p=3の場合のM系列は例えばこう。
0011101
↓ (周期2^3-1=7で同じパターンの繰り返し)
001110100111010011101...
上の桁から3ビットを切り出すと、
001 (1)
_011 (3)
__111 (7)
___110 (6)
____101 (5)
_____010 (2)
______100 (4)

1〜7まで全部出るだろ。これに000だけ追加すればおk。
これだけだと順番がバラバラなので、テーブルと組み合わせる。
(中略)
ビット溢れによるマスクなども組み合わせているが。

そして、TarZさんの見つけた値 0x03F566ED27179461 を使ったのが上のコードです。
http://slashdot.jp/~TarZ/journal/448559

0x03F566ED27179461 =
0000 0011 1111 0101 0110 0110 1110 1101 0010 0111 0001 0111 1001 0100 0110 0001。
以下引用。

このビット列について、6文字を切りだす作業を1文字目から順に行っていくと、以下の通りになる。
2つ目のパターンでsortしてみると、000001から111111まで、すべてのパターンが出現していることが確認できる。
さらにこれに000000を加えればすべて揃うことになる。なんと素晴らしい!

    1    000001
    2    000011
    3    000111
以下略

すごい!すごすぎる!この値はどうやって求めたんだろう?

( y * 0x03F566ED27179461UL ) >> 58 の部分についてもうちょっと書いておくと、y は 2 の n 乗の値、つまり 1, 2, 4, 8,... で掛け算することによって 0x03F566ED27179461 を左シフトすることになります。y が 1 ならシフトなし、2 なら 1bit シフト、4 なら 2bit シフトといったぐあい。これによって桁あふれを起こし、上位の桁を消します(追記。このとき y の値によって何ビット左シフトするか変わるので、y の値によって上位 6bit のビットパターンが変わります)。最後に 58bit (58=64-6) 右シフトして上位 6bit のビットパターンを下位に移動し(下位 bit をマスク)、テーブルを引きます。このとき符号付き整数を使うと符号拡張されてしまい、さらに下位 6bit の AND を取る必要がありますが ulong を使うことでこの手間を回避しています。

参考リンク

英語での解説もみつけました。

おまけ。この問題を解く IEEE754 の float を使ったハックがありますが、それの 64bit版 C# コードも書いてみました。上のコードの10倍ほど遅くなってしまいました><
禁断の unsafe 使えば速いかも?

public static int FloatHack2( long v )
{
    if ( v == 0 ) return 64;

    v = v & -v;

    if ( v == -9223372036854775808 ) return 63;

    float f = ( float ) v;
    var bits = BitConverter.GetBytes( f );
    uint n = BitConverter.ToUInt32( bits, 0 );
    int r = ( int ) ( ( n >> 23 ) - 0x7f );
    return r;
}

(7/6追記) たくさんのアクセスありがとうございます。こんなマニアックな話題なのにまさかのホットエントリー入りに驚いてますw

本当に速いの?ってはてブのコメントに答えると、環境依存してよければ x86 には bsf という機械語があって一命令でカウントできます。x64 の 64bit であれば bsfq。環境依存なしならこのコードは私の知る限り最速です。K8 では bsf よりこっちを使ったほうがよいそうです。詳しくはコメント参照。64回ループするのに比べると手元の環境では10倍程度速いです。速さにこだわっているわけは一つ前の日記あたり。もっと速いのあれば教えてください。上には上がいる!かも?
ちなみに Hacker's delight は機械語にすると何命令実行することになるとか論じてる本ですw 分岐や割り算をできるだけ避けるような、そういう世界。
はてブのわからないってコメントに答えて丁寧に説明してみました。でもわからなくてもいいと思うw

  • 編集履歴
    • 2009/7/6 ( y * 0x03F566ED27179461UL ) >> 58 の説明をちょっと直してみた。リンク1件追加。追記追加。
    • 2009/7/7 「Hacker's delight のコードより4〜5倍速く」だった部分を修正。