原文链接: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为整数的位宽。
2) less-naive solution: 对方法1)进行优化,x&(x-1)可以将最右端的(rightmost)1转化为0
时间复杂度为O(n),n为1的个数。这样在1比较少的时候比第一种方法快很多。
同理,x | (x+1) 可以消除最右端的一个0(转化为1)。因此在0比较多的时候,可以先计算0的个数。下面是wikipedia上给出的版本。跟上面的代码类似,只是把循环展开了。
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的方法。
总体需要 6次shift, 12次and, 6次add 共24次算术运算。
4) 对shift_and_add进行优化
因为具体的分析实在是有点啰嗦,所以先看程序,再来分析。
共需要 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代码如下:
总共需要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方法。还是先看代码:
总共需要 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]:
之前几种方法中的所有操作都可以在寄存器中完成,而查找表法由于无法将查找表放在寄存器中,不可避免地要出现访问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里的内容,就更加向往那个需要很底层优化的时代,想想绞尽脑汁终于得到一个效率更高的解法,简直是太爽了~ 呵呵,不过这个想法估计是实现不了了。
陆续会一点一点把自己的笔记整理出来,希望自己能够坚持下去。