求二进制数中1的个数 (下)

5) 传说中的MIT HAKMEM169方法

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

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方法。还是先看代码：

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
           --------------
            ef gh ij kl
红色标注部分为溢出位，舍掉。橙色标注部分是我们需要的结果，右移24位就可以得到。同理对于64位整数，乘法操作后的结果右移56位就可以得到结果了。

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

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

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) & 0xff];
  return ret;
}

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

三、性能评测

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

naive > less_naive > shift_and_add > shift_add_add_opt > HAKMEM169 > shift_and_add_mul > lut

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

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

参考文献：

[1] Hamming Weight:
[2] Everything上的讨论: %201%20bits
[3] HAKMEM ITEM 169:
[4] Set-bit Population Counts:
[5] Popcount:
[6] HAKMEM 169 and other popcount implementations:
[7] Population count of a 32-bit register: ~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: