素晴らしい。

2009-07-04 - 当面C#と.NETな記録
問題の説明はここまでにして、コードの紹介です。Hacker's delight のコードより4〜5倍速く、そして、イミフ加減が半端じゃない！これ一つで 64bit 値以下のすべての値に対応できます。

でも、実際にどれくらい威力があるか試してみたかったのでCに移植してみた。意外な結果が出ております。

0x03F566ED27179461ULL

まずは黒魔術。より黒魔術っぽくしてみました。

typedef unsigned long long U64;

#define HASH 0x03F566ED27179461ULL
static int ntzhash[64];

void init_ntzhash(){
    U64 h = HASH;
    int i;
    for (i = 0; i < 64; i++, h <<= 1) ntzhash[ h >> 58 ] = i;
}

int ntz_hash(U64 n){
    return n ? ntzhash[ (int)( ( ( n & -n ) * HASH ) >> 58) ] : 64;
}

右シフト

最もナイーブな実装ですね。O(n)

int ntz_shift(U64 n){
    int i;
    for (i = 0; i < 64; i++) if ( (n >> i) & 1 ) return i;
    return 64;
}

二分探索

Hacker's delight P. 85に出てきます。O(log n)、64bitの場合最高8回の条件分岐ということになります。三項演算子を使った方が見やすいのでそうしました。

int ntz_bsrch( U64 n ) {
    n &= -n;
    return n == 0 ? 64
    : n == 0x80000000ULL ? 31
    : (n <  0x80000000ULL ? n == 0x8000ULL ? 15
    /* 60行ほど中略 */
               : (n <  0x4000000000000000ULL ? 61
              : 63))))));
}

こんなひちめんどくさいものは手で書いていられないので、以下のPerl Scriptに自動生成してもらいました(要use64bitint)。

sub foo{
    my ($head, $tail) = @_;
    my $mid = int(($head + $tail) / 2);
    return $mid if $head >= $tail;
    my $ret = sprintf "n == 0x%xULL ? %d\n", 2**$mid, $mid;
    $ret .= 
        sprintf ": (n <  0x%xULL ? %s\n: %s)", 2**$mid,
        foo($head, $mid-1),    foo($mid+1, $tail);
    $ret =~ s/\n+/\n/g;
    $ret;
}

print foo(0,63)

frexp

n & -nで、最低位ビットだけ立てられるということは、二のベキ乗にすでになっているわけです。そこからどうやってビット数にマップするかが、黒魔術のすごいところなのですが、もう少しけれんみのないやり方として、frexp()を使うというやり方もあります。

#include <math.h>

int ntz_frexp(U64 n){
    int z;
    if ( n == 0 ) return 64;
    frexp( (double)( n & -n ), &z);
    return z - 1;
}

x86のbsf命令

で、x86のbsf命令を使うやり方です。どのx86でも動くよう、ここでは32bit以上なのか否かで条件分岐させています。

int ntz_bsf(U64 n){
    if (n == 0) return 64;
    unsigned int lo = n & 0xffffffff;
    int z;
    if (lo){
        __asm__("bsf %1, %0;" :"=r"(z) :"r"(lo));
    return z;
    }else{
        n >>= 32;
        __asm__("bsf %1, %0;" :"=r"(z) :"r"(n));
        return z + 32;
    }
}

まとめてベンチマーク

で、この５つをまとめてベンチマークしてみました。引数に繰り返し回数を指定するとベンチマーク、何もなければテストになります。Codepadでも動きました。

http://codepad.org/9UXR2cOu 
#include <sys/time.h>

double time_hires(){
    struct timeval now;
    gettimeofday(&now, NULL);
    return now.tv_sec + now.tv_usec/1e6;
}

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

#define D63 ((U64)1 << 63)

typedef struct { char *name; int (*func)(U64); } nfpair;

static nfpair ntz[] = {
    {"ntz_shift", ntz_shift },
    {"ntz_bsrch", ntz_bsrch },
    {"ntz_frexp", ntz_frexp },
    {"ntz_hash",  ntz_hash  },
    {"ntz_bsf",   ntz_bsf   }
};

int main (int argc, char **argv){
    U64 n;
    int i, j, count = 0;
    double timer;
    if (ntzhash[63] != 6) init_ntzhash();
    if (argc > 1){ /* benchmark */
        count = (int)strtod(argv[1], NULL);
        for (j = 0; j < 5; j++){
            timer = time_hires();
            for (i = 0; i < count; i++){
                assert(ntz[j].func(D63) == 63);
            }
            printf("%9s: %gs\n", ntz[j].name, time_hires() - timer);
        }
    }else{ /* test */
	printf("%d\n", ntz_bsrch(D63));
        for (j = 0; j < 5; j++){
            for (i = 0, n = 1 ; i <= 64; i++, n += n){
                assert(ntz[j].func(n) == i);
                printf("ok %d # %s(0x%016llx) = %d\n", 
                       ++count, ntz[j].name, n, ntz[j].func(n));
            }
        }
        printf("1..%d\n", count);
    }
    return 0;
}

で、以下が結果です。

Core 2 Duo 2.4GHz, Mac OS X v10.5.7 
% gcc -W -Wall -O6 ntz.c && ./a.out 1e8
ntz_shift: 16.1367s
ntz_bsrch: 1.8118s
ntz_frexp: 2.79462s
 ntz_hash: 1.18244s
  ntz_bsf: 0.475164s
Dual Xeon 2.66GHz, FreeBSD 7.2 
% gcc -W -Wall -O6 ntz.c && ./a.out 1e8
ntz_shift: 25.2944s
ntz_bsrch: 1.65888s
ntz_frexp: 7.1559s
 ntz_hash: 1.50684s
  ntz_bsf: 0.905658s

ベンチマークで使っている数値は、O(c)でないアルゴリズムで最悪になるよう、1<<63にしてあります。見ての通り、二分探索と黒魔術の差が以外と出てません。64bitの乗算のコストかなあ....

おまけ

Perlの一行野郎だとこんなでしょうか。0の場合も対応しています。一行野郎の割には素直すぎる....

sub ntz{ length(sprintf("%064b", shift) =~ /(0+)$/ && $1) }

Dan the Bit Counter