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

一、问题描述：求一个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解法非常精彩，让人看完以后不禁拍案叫绝又扼腕叹息，为什么我就想不到这种方法呢 T.T。

二、解法分析：

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

typedef unsigned __int64 uint64
int popcount_naive(uint64 x)
{
    int count = 0;
    for(; x; x>>1)
        ++count;
    return count;
}

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

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上给出的版本。跟上面的代码类似，只是把循环展开了。

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的方法。

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进行优化

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

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就可以了。