Pop Count Problem (二进制数中1的个数)

原文链接:http://blog.ibread.net/375/pop-count-problem/

以下文章中有大红色背景的文字是我转载时做的注释。


一、问题描述:求一个N位整数x的二进制表示中1的个数,越快越好。

据说这是道很著名的面试题。原题是说如何在常量时间内算出32位整数的二进制表示中1的个数。实际上这么问是有漏洞的,因为按照最笨的方法,一个一个数,也不过32次,当然是常量时间。但如果你这么告诉面试官,十有八九会被骂白痴。因为大家一般都会认为32次实际上是O(N)的做法。所以这么问可能会更好一些,求一个N位整数x的二进制表示中1的个数,越快越好。

实际上这个问题叫做Hamming weight[1],或者叫做population count以及pop count。看到Hamming大家基本上可以恍然大悟一下了,这玩意可以用来计算海明距离。另外在信息论、编码学等也有很多应用。

对于这个问题的解法,Google一下popcount 或者count bits 1会有一大把。基本上无非是三大类,直观的挨个数O(N),分治法O(lgN),查表法O(1)。但实际上对于这种问题,单单考虑时间复杂度意义不大。尤其是在问题规模较小的情况下,比如说32位整数或者64位整数。比方说,对于查表法,如果查找表太大导致程序运行时不得不访存,这个开销可能要比前两种方法都要大。下面分别总结一下,以64位整数为例,源程序和解析都给出来了。这个题目虽然简单,但很多解法,尤其是HAKMEM解法非常精彩,让人看完以后不禁拍案叫绝又扼腕叹息啊。

二、解法分析:

1) naive solution:依次移位,数1的个数,时间复杂度为O(N)。这里的N为整数的位宽。

1
2
3
4
5
6
7
8
typedef unsigned __int64 uint64;
int popcount_naive(uint64 x)
{
    int count = 0;
    for(; x; x>>=1) // 原文是for(; x; x>>1),怀疑少了一个'='号
++count; return count;}

2) less-naive solution: 对方法1)进行优化,x&(x-1)可以将最右端的(rightmost)1转化为0

1
2
3
4
5
6
7
int pop_count_lnaive(uint64 x)
{
    int count = 0;
    for (; x; ++count)
        x&=x-1;
    return count;
}

时间复杂度为O(n),n为1的个数。这样在1比较少的时候比第一种方法快很多。
同理,x | (x+1) 可以消除最右端的一个0(转化为1)。因此在0比较多的时候,可以先计算0的个数。下面是wikipedia上给出的版本。跟上面的代码类似,只是把循环展开了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
const uint64 hff = 0xffffffffffffffff;
#define f(y) if ((x &= x-1) == 0) return y;
int popcount_5(uint64 x) {
    if (x == 0) return 0;
    f( 1) f( 2) f( 3) f( 4) f( 5) f( 6) f( 7) f( 8)
        f( 9) f(10) f(11) f(12) f(13) f(14) f(15) f(16)
        f(17) f(18) f(19) f(20) f(21) f(22) f(23) f(24)
        f(25) f(26) f(27) f(28) f(29) f(30) f(31) f(32)
        f(33) f(34) f(35) f(36) f(37) f(38) f(39) f(40)
        f(41) f(42) f(43) f(44) f(45) f(46) f(47) f(48)
        f(49) f(50) f(51) f(52) f(53) f(54) f(55) f(56)
        f(57) f(58) f(59) f(60) f(61) f(62) f(63)
        return 64;
}
//Use this instead if most bits in x are 1 instead of 0
#define f(y) if ((x |= x+1) == hff) return 64-y;

3) shift_and_add 时间复杂度O(logN)
这个方法主要利用两点,一是根据最直观的做法,把每一位相加后就可以得到结果,这个过程我们可以利用分治+并行加法来优化;二是对于n位整数,最多有n个1,而n必定能由n位二进制数来表示,因此我们在求出某k位中1的个数后,可以将结果直接存储在这k位中,不需要额外的空间。
以4位整数abcd为例,最终结果是a+b+c+d,循环的话需要4步加法。
1′ 每一位一组[a][b][c][d],相邻两组相加得到[a+b][c+d],结合成两位一组。所有的加法可以像这样并行进行:

[0] [b] [0] [d]
[0] [a] [0] [c]
---------------
[e   f] [g   h]

其中ef=a+b, gh=c+d。而0b0d = (abcd) & 0101, 0a0c = (abcd)>>1 & 0101
2′ 两位一组: [ef][gh],相邻两组相加,结合成四位一组。

[0 0] [g h]
      [e f]
-----------
  [i j k l]

其中ijkl = ef + gh,00gh = (efgh) & 0011, 00ef = (abcd)>>2 && 0011。这样得出的ijkl就是最终结果。这样只需要log(N)步。
按照这个方法类推,我们很容易可以得出计算uint64的方法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
typedef unsigned __int64 uint64; //assume this gives 64-bits
const uint64 m1 = 0x5555555555555555; //binary: 0101...
const uint64 m2 = 0x3333333333333333; //binary: 00110011..
const uint64 m4 = 0x0f0f0f0f0f0f0f0f; //binary: 4 zeros, 4 ones ...
const uint64 m8 = 0x00ff00ff00ff00ff; //binary: 8 zeros, 8 ones ...
const uint64 m16 = 0x0000ffff0000ffff; //binary: 16 zeros, 16 ones ...
const uint64 m32 = 0x00000000ffffffff; //binary: 32 zeros, 32 ones
 
int popcount_1(uint64 x) {
 x = (x & m1 ) + ((x >> 1) & m1 ); 
 x = (x & m2 ) + ((x >> 2) & m2 ); 
 x = (x & m4 ) + ((x >> 4) & m4 ); 
 x = (x & m8 ) + ((x >> 8) & m8 ); 
 x = (x & m16) + ((x >> 16) & m16); 
 x = (x & m32) + ((x >> 32) & m32); 
 return x;
}

总体需要 6次shift, 12次and, 6次add 共24次算术运算。

4) 对shift_and_add进行优化

因为具体的分析实在是有点啰嗦,所以先看程序,再来分析

1
2
3
4
5
6
7
8
9
10
int popcount_2(uint64 x)
{
    x -= (x >> 1) & m1;
    x = (x & m2) + ((x >> 2) & m2);
    x = (x + (x >> 4)) & m4;
    x += x >> 8;
    x += x >> 16;
    x += x >> 32;
    return x & 0x7f;
}

共需要 6次shift, 5次and, 5次add, 1次sub,共17次算数运算。下面具体来看这7个与运算是如何省掉的。

1′ 第一步要做的是将2bits的整数ab (=2*a+b)变成 a+b。不难看出,(2*a+b) – a = a+b。也就是说: x – (x>>1) & m1就是我们要的结果。这样省掉了1步与运算

2′ 先与后加=>先加后与,再省掉1步与运算。

之所以先与后加,是为了去掉不必要的干扰。比如说对于上面所讲的例子:

[a] [b] [c] [d]
[0] [a] [b] [c]
—————
[e   f] [g   h]

我们所关心的结果只是奇数组(d+c)和(b+a)的结果,这里我们称之为有效组; 红色标注的部分(c+b)和(a+0)的结果是没有用的(上例中直接置为0),我们称之为无效组。每一步的操作,实际上是将有效组相加,结果存放在有效组和无效组连接所组成的新组中

如果不通过and操作将abcd=>0b0d, 0abc=>0a0c,那么在做加法的时候无效组会出来捣乱。第一组的结果(d+c)如果有进位,就会被第二组的结果(c+b)所干扰。同理,第二组(c+b)的进位也可能会干扰到第三组的结果(b+a)。这样的无效组我们称之为干扰组,因此在加法操作进行之前有必要利用掩码将干扰组屏蔽

通过上面的分析我们可以知道,干扰组的出现是由于两个k位分组的最大和(2k)超出了k位二进制数能表示的范围,从而不可避免地产生进位,干扰相邻组的结果那么如果不会产生进位呢,就没有必要在加法操作之前利用掩码将干扰位过滤掉了,过滤操作可以延迟到加法操作之后进行。也就是说,我们可以先加后与,省掉一个与操作。(就为了省掉一个与运算,还真是折腾。。。)

不难看出,从4位分组开始,两个分组的最大和8就可以仅用原来的4位表示了。因此m4, m8, m16, m32都可以延迟到加法操作之后。

3′ 只加不与,再省掉1步与运算

之所以加法操作之后还是要进行与运算,是为了防止将脏数据带入下一步计算中。比如说:

[a1 a2 a3 a4] [b1 b2 b3 b4] [c1 c2 c3 c4] [d1 d2 d3 d4]
[a1 a2 a3 a4] [b1 b2 b3 b4] [c1 c2 c3 c4]
——————————————————-
[e1 e2 e3 e4] [f1 f2 f3 f4] [g1 g2 g3 g4] [h1 h2 h3 h4]

很明显,红色标注的部分的无效组并不会产生进位,从而不需要在加法操作之前过滤。那么如果加法操作之后还是不过滤呢?如果继续将这部分脏数据带到下一步运算中呢?

[e1 e2 e3 e4 f1 f2 f3 f4] [g1 g2 g3 g4 h1 h2 h3 h4]
[e1 e2 e3 e4 f1 f2 f3 f4]
—————————————————-

[h1 h2 h3 h4] + [f1 f2 f3 f4]最大是16, 很显然4位已经存不下了,但上一步操作中的脏数据。也就是说,无效组即便在当前步骤中不会捣乱,也可能在下一步中出来捣乱,我们可以称之为延迟干扰。因此及时斩草除根免除后患是必要的。

那么如果无效组在今后所有的步骤中都不会形成延迟干扰呢?也就是说,之后每一步的中间结果都只用k位就能存下,不再需要进位。这样一来,每次计算只需要关心分组中后k位的结果就可以了,无效位即使不清除,也不会造成干扰,因此加法操作后的与运算也可以省掉

不难看出,64位整数最多有64个1,而64只需要7位就可以存下。也就是说,从8位一组两两相加开始,我们只需要关心每组中后7位的运算结果,而且这后7位的结果不会受到任何干扰。因此m8, m16, m32全部都可以不要,只需要在最后return x& 0x7f就可以了。

5) 传说中的MIT HAKMEM169方法

HAKMEM 是1972年由MIT人工智能实验室的一帮牛人(当然也包括RMS)发布的一本备忘录,里面都是些很无敌的算法,主要的目标是更快更好地进行数学运算。其中的第169个算法,就跟popcount有关,原始的代码[3]是用汇编写的,翻译成C代码如下:

1
2
3
4
5
6
7
8
9
10
typedef unsigned __int32 uint32;
const uint32 m1=033333333333;
const uint32 m2=011111111111;
const uint32 m3=030707070707;
 
int popcount_hakmem(uint32 x) {
    x -= (x>>1 & m1 + x>>2 & m2);
    x += x>>3 & m3;
    return x % 63;
}

总共需要3次shift, 3次and, 2次add, 1次mul, 1次mod 共10次算数运算。这是32位整数的版本,改成适用于64位整数的版本也很简单。主要思想也是分治以及并行加法。只是与原始的shift_and_add有两点不同:

一是第一步的时候是每位一组,相邻三组相加。比如说对于x=(abc)2=4*a+2*b+c,x>>1&m1= 2*a+c,x>>2&m2=a,于是 x – x>>1&m1 – x>>2&m2 = a+b+c。
二是在第二步做完之后,这时是每6位自成一组,以12位整数为例,x=(ab)64 = a*64 + b。很明显, a*64+b 
≡ a+b mod 63。也就是说x与a+b关于模63同余。同理,对于64位整数我们也可以这么处理。

6) 再次优化shift_and_add方法,用乘法操作代替加法操作

利用CPU中的基本运算部件乘法器,我们可以继续优化shift_and_add方法。还是先看代码:

1
2
3
4
5
6
7
const uint64 h01 = 0x0101010101010101;
int popcount_3(uint64 x) {
    x -= (x >> 1) & m1;
    x = (x & m2) + ((x >> 2) & m2);
    x = (x + (x >> 4)) & m4;
    return (x * h01)>>56;
}

总共需要 4次shift, 4次and, 2次加法,1次减法,1次乘法共12次算数运算。

可以看出,第三步以后的操作全部都被省略了。这是因为从这步开始,我们需要的操作是8位一组,8个组的数依次相加,就可以得到最终结果。这恰恰可以通过通过乘法器来做到。以32位整数为例,共有四组:

            0a 0b 0c 0d
        *   01 01 01 01
           -------------
            0a 0b 0c 0d
         0a 0b 0c 0d
      0a 0b 0c 0d
   0a 0b 0c 0d
           --------------
            0a gh ij kl

红色标注部分为溢出位,舍掉。橙色标注部分是我们需要的结果右移24位就可以得到。同理对于64位整数,乘法操作后的结果右移56位就可以得到结果了

注意如果CPU执行乘法操作指令比较慢的话,这样优化可能会适得其反。但一般来说不会这样,比如说AMD在两个时钟周期里就可以完成乘法运算。

6) 最笨最快的查表法:
很多情况下我们都可以以牺牲空间的代价来得到时间上的最优解。比如说对于这道题目,我们枚举8位整数的二进制中1的个数,组成一个查找表,看代码[4]:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
const int lut[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
};
 
static inline int popcount_lut(uint64 x) {
    int i, ret = 0;
 
    for (i = 0; i < 64; i += 8)
        ret +=lut[(x >> i) &amp; 0xff];
    return ret;
}

之前几种方法中的所有操作都可以在寄存器中完成,而查找表法由于无法将查找表放在寄存器中,不可避免地要出现访问cache或内存的操作。但由于我们采用一个循环频繁访问查找表,在所有的实现中查找表都应该被锁定在cache中,访问开销很小[4]。因此从后面的评测我们可以看到,这种方法的效率最高。

三、性能评测

实际上我在总结这几种方法的时候,基本上是按照性能逐渐提升的顺序来的。也就是说对于时间开销:

对于很酷的HAKMEM169,由于用到了耗时的mod操作,在大多数平台上都不如shift_and_add的乘法版本速度快,当然就更比不上我们最笨最好的查找表大法了。KISS准则再次被证明是有效的,哦也。具体的评测数据,可以参考Bart Massey的结果[5]以及感言,或者看这个更加全面详尽啰嗦的评测[6]。

实际上,Intel已经将这个数1操作封装成一个指令popcnt,在Itanium中就已提供。既然CPU提供原生支持,那想必应该是最快的了(据说1个时钟周期就可以[10])。

参考文献:

[1] Hamming Weight: http://en.wikipedia.org/wiki/Hamming_weight
[2] Everything上的讨论: http://everything2.com/index.pl?node=Counting%201%20bits
[3] HAKMEM ITEM 169: http://www.inwap.com/pdp10/hbaker/hakmem/hacks.html#item169
[4] Set-bit Population Counts: http://www.technovelty.org/code/arch/popcount.html
[5] Popcount: http://wiki.cs.pdx.edu/forge/popcount.html
[6] HAKMEM 169 and other popcount implementations: http://www.dalkescientific.com/writings/diary/archive/2008/07/03/hakmem_and_other_popcounts.html
[7] Population count of a 32-bit register: http://www.df.lth.se/~john_e/gems/gem002d.html
[8] 《编程之美》读书笔记——“求二进制数中1的个数: http://bvcat007.javaeye.com/blog/203577
[9] Henry S. Warren, Jr. Chapter 10. The Quest for an Accelerated Population Count. Beautiful Code.
[10]GMP: http://gmplib.org/

后记:

这篇本来是09年寒假的时候发在live space上的,原文链接在这里。时间过得太快,现在live space都已经不存在了,我把它迁移到了wordpress上。已经记不太清楚当时总结这些的心情了,只隐约记得当时应该是在成都,很冷。心情应该也不会好到哪里去,因为在纠结是否出国。最开始知道这个题目是LSH师兄去Adobe笔试碰到的题目,回来以后大家讨论了解法,便一直想找时间总结下来。我做事情有个毛病,就是尽可能地要尽善尽美,所以就一直拖一直拖,直到寒假有了整块的时间才写下来,而且很罗嗦。

隐约还记得第一次看到HAKMEM解法时的惊讶,那时就想,原来还可以这么做。后来又陆续看了一些Hakers’ delight里的内容,就更加向往那个需要很底层优化的时代,想想绞尽脑汁终于得到一个效率更高的解法,简直是太爽了~ 呵呵,不过这个想法估计是实现不了了。

陆续会一点一点把自己的笔记整理出来,希望自己能够坚持下去。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值