原帖
http://topic.youkuaiyun.com/u/20110426/15/889BF7F1-409B-4F15-91DF-1E21C5A14743.html
http://topic.youkuaiyun.com/u/20110427/13/6F8371C6-897D-44A3-9DA7-E6FD3B4DE761.html
问题:40位位计数.
设计一个算法,输入一个64位整数,求其低40位中1的个数.已知高24位必定为0.
算法1:循环.
int bitcount_loop(UINT64 s)
{
int c;
for (c = 0; s; s >>= 1)
{
c += s & 1;
}
return c;
}
众所周知,这种算法很低效的.
算法2:依次取出尾部的1.
int bitcount_trail(UINT64 s)
{
int c;
for (c = 0; s; s &= s - 1)
{
c++;
}
return c;
}
比算法1循环次数少了点,但仍不够快.
s &= s - 1;意为把s最低的1位清0.
算法3:查表.
int bitcount_table(UINT64 s)
{
static const BYTE t[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};
BYTE* n = (BYTE*)&s;
return t[n[0]] + t[n[1]] + t[n[2]] + t[n[3]] + t[n[4]];
}
预先计算出8位的位计数,把s的40位分成5个字节分别查表在相加.
查表意味着空间换时间,在cpu不够先进的年代,这种算法是最快的.
但现代的cpu做数学运算很快的,相比之下访问内存慢.
每16查表也可以,但这样需要一个64k的大表,使高速缓存命中率降低,效率不一定能提高.
32位的算法.
int bitcount32(UINT x)
{
x -= x >> 1 & 0x55555555;
x = (x >> 2 & 0x33333333) + (x & 0x33333333);
x = (x >> 4) + x & 0x0F0F0F0F;
x += x >> 8;
x += x >> 16;
return x & 63;
}
这是32位位计数的经典算法,是并行算法的一个优化版本.在大部分机器上,这种算法是最快的.
如果不明白,请看<<高效程序的奥秘>>.
算法4:调用2次32位算法.
int bitcount_use32(UINT64 s)
{
return bitcount32(UINT(s)) + bitcount(UINT(s>>32));
}
分别对s的低32位和高32位使用bitcount32.
其中UINT(s>>32)是s的高32位.一个合格的编译器不会真的执行移位运算,而是直接取高位.
我个人的习惯:如果一个表达式可在编译时求值,则其中的运算符前后不加空格.
算法5:
int bitcount_insert(UINT64 s)
{
UINT c = UINT(s), d = UINT(s>>32);
c -= c >> 1 & 0x55555555;
c += d + (d << 15) & 0x55555555; // 强力插入
c = (c >> 2 & 0x33333333) + (c & 0x33333333);
c = (c >> 4) + c & 0x0F0F0F0F;
c += c >> 8;
c += c >> 16;
return c & 63;
}
前几种算法未利用"高24位必定为0"这一条件.
这种算法除了标有"强力插入"这一行,跟32位算法完全相同.
为什么d + (d << 15) & 0x55555555,我用加不用或?因为在vc6中,用或导致多一条与指令.
此算法对最多47位有效.如果用或,最多48位有效.
我的程序中只用到40位,我就用加了.
此行作用是将d的0,2,4,6位保持不变,将d的1,3,5,7位移动到16,18,20,22位中.
上一行执行完后,c中的每两位看做一个单元,其值可能是00,01,10.
给每个单元加上一个单独的位,不会溢出.
此后的计算,跟32位完全相同.
_________________________
最后一段我说的不清楚。
d + (d << 15) & 0x55555555
vc6会如实生成机器码。
(d | (d << 15)) & 0x55555555
vc6会“优化”为(d & 0x55555555) | ((d & 0xaaaa) << 15)
速度稍有变慢。我很奇怪编译器这样优化道理何在。
=========================================================
设计一个算法,输入一个64位整数s,这个整数可表示为s = 1 << c; 其中c < 40. 求c的值.
即:求s中唯一的1的位置.
算法1:
int bitposition_loop(UINT64 s)
{
for (int c = 0; c < 40; c++)
if (s == 1ULL << c)
return c;
}
最直观的,也是最慢的算法.
算法2:
int bitposition_count(UINT64 s)
{
return bitcount(s - 1);
}
s-1可生成c个连续的1.对其计数即可.
算法3:
int bitposition_rem(UINT64 s)
{
static const BYTE T[53] = { // 余数-指数表, c == (1 << t[c]) % 53
-1, 0, 1, 17, 2, 47, 18, 14, 3, 34, 48, 6, 19, 24, 15, 12, 4, 10,
35, 37, 49, 31, 7, 39, 20, 42, 25, 51, 16, 46, 13, 33, 5, 23, 11, 9,
36, 30, 38, 41, 50, 45, 32, 22, 8, 29, 40, 44, 21, 28, 43, 27, 26};
return T[(UINT(s>>32) * 42 + UINT(s)) % 53];
}
每一个(1 << c) % 53都有一个独特的余数.知道余数就能查表得出c.
除数53不是随便选的,跟c的范围有关.比如c < 8, 除数要选11. c < 32, 除数要选37.
为了避免64位除法,使用了一些技巧.设s的高32位为a, 低32位为b,
s = a * 0x100000000 + b;
而0x100000000 = 4294967296 = 81037118 * 53 + 42;
s % 53 = (a * 0x100000000 + b) % 53
= (a * (81037118 * 53 + 42) + b) % 53;
= (a * 81037118 * 53 + a * 42 + b) % 53
= (a * 81037118) * 53 % 53 + (a * 42 + b) % 53
= (a * 42 + b) % 53;
尽管如此,除法还是一个相当慢的操作.所以,此算法效率不高.
算法4:
int bitposition_float(UINT64 s)
{
float f = (float)(INT64)s;
return ((int&)f >> 23) - 127;
}
算法5:
int bitposition_double(UINT64 s)
{
double d = (double)(INT64)s;
return UINT((UINT64&)d >> 52) - 1023;
}
这两个算法是把64位整数转化为浮点数,然后取阶码.
理解这个算法要对浮点数格式有所了解.
IEEE754规定,浮点数由符号,阶码,尾数3部分组成.
float是32位浮点数,其中符号占1位,阶码占7位,尾数占23位.
double是64位浮点数,其中符号占1位,阶码占11位,尾数占52位.
long double是80位浮点数,什么占多少位忘记了.
阶码是用移码表示的.float的移位是127, double是1023.
还有一个奇怪的语法(int&)f,他表示将f所在的内存中的"东西"看做int,相当于*(int*)&f.
(UINT64&)d同理.
(double)(INT64)s看似奇怪,但没有(INT64)不行,编译器会说不能把UINT64转化为float.
理由是long double放不下64位有效数字.我不能理解.
这两个算法生成的机器码包含以下几个指令:fild,fstp,shr,sub.指令虽少,前两个指令是比较慢的.
算法6:
int bitposition_bin(UINT64 s) //23834
{
int c, a;
if (UINT(s))
{
c = 0;
a = UINT(s);
}
else
{
c = 32;
a = UINT(s>>32);
}
if (a & 0xffff0000) c += 16;
if (a & 0xff00ff00) c += 8;
if (a & 0xf0f0f0f0) c += 4;
if (a & 0xcccccccc) c += 2;
if (a & 0xaaaaaaaa) c += 1;
return c;
}
2分查找法.指令较多,但都是快速指令.速度还凑合.
算法7:
int bitposition_mul(UINT64 s)
{
static const BYTE T[32] = {
0, 1, 2, 6, 3, 16, 7, 21, 4, 19, 17, 11, 13, 8, 22, 26,
31, 5, 15, 20, 18, 10, 12, 25, 30, 14, 9, 24, 29, 23, 28, 27};
if ((UINT)s)
return T[(UINT)s * 0x046B29DF >> 27];
else
return T[UINT(s>>32) * 0x046B29DF >> 27] + 32;
}
0x046B29DF,这个神奇的数,他的2进制可表示为00000100011010110010100111011111
仔细观察可发现,这个序列中包含了从00000到11111中所有可能的组合.这样的序列称为DeBruijn序列.
0x046B29DF * s >> 27 = 0x046B29DF * (1 << c) >> 27 = 0x046B29DF >> 27 - c;
就像一张写字的纸条,透过一个小孔看纸条,根据看到的文字就能判断纸条的位置.
顺便提一句,DeBruijn序列和Gray码在位置传感器上都有应用.
这种算法仅需一次乘,一次移位,一次查表,在乘法比较快的现代计算机上速度还是不错的.
算法8:
int bitposition_bsf(UINT64 s)
{
if ((UINT)s)
{
__asm bsf eax, s
}
else
{
__asm bsf ecx, [s+4]
__asm lea eax, [ecx+32]
}
}
bsf,一个被遗忘的指令.
格式:bsf dst, src
功能:查找src最低的"1"位,其位置保存在dst中,影响标志位ZF.其中src必须是寄存器.
多么美妙的指令,正是我想要的!
且慢!bsf,由于其太过设备相关,几乎没有高级语言直接支持它.
而得不到支持的指令,硬件设计师肯定不会花时间优化它.得不到优化就更得不到支持...
唉,我这忧国忧民的习惯啥时能改改呢?在bsf从指令集消失之前,看一下他的表现吧!
结果怎样?不提也罢...
相比bitcount,bitposition的几个算法很难明显的分出优劣.一些算法可能在某些机器上胜出,
而在另一些机器上表现不佳.但如果本来就要查表的话,比如
sometype somefunction(UINT64 s)
{
static const sometype F[40] = {......};
return F[bitposition(s)];
}
此时bitposition_mul可以无缝的融入somefunction中.
sometype somefunction(UINT64 s)
{
static const sometype F[40] = {......}; // 将F中的元素按DeBruijn序列重新排序
if ((UINT)s)
return F[(UINT)s * 0x046B29DF >> 27];
else
return (F+32)[UINT(s>>32) * 0x17000000 >> 29];
}
其中0x17000000是一个3阶DeBruijn序列,他能根据3个bit确定移位的数目,限定移位数<8.
省去了一次查表,速度上还是略微占优的.当然这样程序更加晦涩了.