Hamming Weight的算法分析

看代码时遇到一个求32bit二进制数中1的个数的问题，感觉算法很奇妙，特记录学习心得于此，备忘。

计算一个64bit二进制数中1的个数。

解决这个问题的算法不难，很自然就可以想到，但是要给出问题的最优解，却很有难度。

通常，最容易想到的算法是除余法，继而考虑到除法的代价较高，而且除数是2，会想到使用向右移位来代替除法，并使用&0x1操作来取末位的值，这样提高了算法的效率。然而，这样仍然进行了63次&操作、63次移位操作和63次+操作。若假设字长大小不限，记作N，那么上述算法的时间复杂度都为O(N)。

当然，还有更优的算法。

这个问题其实是HammingWeight的一个应用，又叫做populationcount,popcountorsidewayssum。HammingWeight详见http://en.wikipedia.org/wiki/Hamming_weight，以下部分内容取自维基百科。

Hammingcode是指一个字串中非0符号的个数（TheHamming weight of a stringis the number of symbols that are different from the zero-symbol ofthealphabetused.）。应用到2进制符号序列中来，即二进制串中1的个数就是该串的Hammingcode.那么上述的问题即转换成求解字串的Hammingcode的问题。

下面对维基百科上给出的算法，进行分析。算法使用c语言实现。

//types and constants used in the functions below
 
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
const uint64 hff = 0xffffffffffffffff; //binary: all ones
const uint64 h01 = 0x0101010101010101; //the sum of 256 to the power of 0,1,2,3...
 
//This is a naive implementation, shown for comparison,
//and to help in understanding the better functions.
//It uses 24 arithmetic operations (shift, add, and).
int popcount_1(uint64 x) {
    x = (x & m1 ) + ((x >>  1) & m1 ); //put count of each  2 bits into those  2 bits 
    x = (x & m2 ) + ((x >>  2) & m2 ); //put count of each  4 bits into those  4 bits 
    x = (x & m4 ) + ((x >>  4) & m4 ); //put count of each  8 bits into those  8 bits 
    x = (x & m8 ) + ((x >>  8) & m8 ); //put count of each 16 bits into those 16 bits 
    x = (x & m16) + ((x >> 16) & m16); //put count of each 32 bits into those 32 bits 
    x = (x & m32) + ((x >> 32) & m32); //put count of each 64 bits into those 64 bits 
    return x;
}

分析：popcount1是下面算法的基础，理解了这个算法的思想，下面的算法不过就是此算法的局部优化罢了。
首先，理解这样一个事实，64bit的二进制串中最多有64个1，而0~64内的值必然可以使用该串的低8位来表示（2^8>64）。2^2>2,那么2bit的串中的1的个数必然可以用这两位来表示。
我们先简化成8bit的串，来描述算法的基本思想。
使用abcdefgh来代表一个8bit的2进制串，其中a,b,c,d,e,f,g,h属于集合{0,1}
那么算法求解的目标输出是out=a+b+c+d+e+f+g+h
对应到上面代码中的第一步来说，x = (x & m1 ) + ((x >> 2) & m1 )，
x&m1 = 0b0d0f0h
(x>>2)&m1 = 0a0c0e0g
求和得到：[a+b]2[c+d]2[e+f]2[g+h]2,这里[x]2 表示2位的二进制，其值=x(x表示10进制的值)。如果对应到64bit的串，那么这里将有32个2-bit的组合，即将64bit两两一组，并使用其来表示自身包含的1的个数。

代码的第二步：x = (x & m2 ) + ((x >> 4) & m2 ),同样使用8bit串来简化描述。
x&m2 = 00[c+d]200[g+h]2
(x>>4)&m2 = 00[a+b]2 00[e+f]2
求和得到：[a+b+c+d]4[e+f+g+h]4

第三步： x = (x & m4 ) + ((x >> 4) & m4 );
x&m4 = 0000[e+f+g+h]4
(x>>4)&m2 = 0000[a+b+c+d]4
求和得到：[a+b+c+d+e+f+g+h]8
至此问题得解。对于64bit的串，则如代码所示还要多进行3步。

到这里可以很清楚的看到，算法是使用了分治的思想，每步将问题划分成子问题，然后合并来减小问题的规模，求解问题的过程像是一棵倒置的二叉树。先将n位的二进制相邻的两位两两分为一组，并巧妙的利用移位和掩码来使其利用自身来表示所得到的和，这样从宏观上来看，问题就被简化成规模为n/2bit(这里的bit其实已经是虚指了，其实理解为unit更好)的问题求解了，同样的，继续两两划分成一组分治求解。经过lg（n）步，得到最终的解。
由以上分析可见，算法的复杂度为O（lgn）。
对于64位的字串来说 ，只使用了24次算数操作，比起前面的算法来说要明显减少了。

//This uses fewer arithmetic operations than any other known  
//implementation on machines with slow multiplication.
//It uses 17 arithmetic operations.
int popcount_2(uint64 x) {
    x -= (x >> 1) & m1;             //put count of each 2 bits into those 2 bits
    x = (x & m2) + ((x >> 2) & m2); //put count of each 4 bits into those 4 bits 
    x = (x + (x >> 4)) & m4;        //put count of each 8 bits into those 8 bits 
    x += x >>  8;  //put count of each 16 bits into their lowest 8 bits
    x += x >> 16;  //put count of each 32 bits into their lowest 8 bits
    x += x >> 32;  //put count of each 64 bits into their lowest 8 bits
    return x & 0x7f;
}

popcount2在popcount1的基础上进行了优化。
第一步基于了这样一个事实:ab-0a得到的值为ab中1的个数。
简单证明：若a为0，那么0a=0,减0无变化，那么b就是结果。
若a位1，那么只有两种情况，10-01 = 01, 11-01 = 10.都符合上述事实。
这样x -= (x >> 1) & m1和 x = (x & m1 ) + ((x >> 1) & m1 )的结果相同，却节省了1个操作。（这里我有个疑问，有符号数使用补码进行减法操作等于加法操作，效率相当，然而这里x为无符号数，即原码加减法，原码的减法在机器级如何实现，即一个源码减法的操作的代价与加法和与操作的代价和比较，哪个更大？有时间的话要去看下原码减法的实现）

第二步第三步同popcount1，此时x=[a]8[b]8[c]8[d]8[e]8[f]8[g]8[h]8

第四步后x = [H8|a+b]16[H8|c+d]16[H8|e+f]16[H8|g+h]16,这里H8代表高8位，由于我们不关心高8位的值（当然H的值是明显知道的），这里就用H代替。由于使用低8位完全可以表示0～64范围内的值，因此不用担心低八位溢出。

同理，第五步后x=[H24|a+b+c+c]32[H24|e+f+g+h]32
第六步后x=[H56|e+f+g+h]64
第七步使用掩码0x7f获得低8位的值(0xff效果应该一样的吧..?)

//This uses fewer arithmetic operations than any other known  
//implementation on machines with fast multiplication.
//It uses 12 arithmetic operations, one of which is a multiply.
int popcount_3(uint64 x) {
    x -= (x >> 1) & m1;             //put count of each 2 bits into those 2 bits
    x = (x & m2) + ((x >> 2) & m2); //put count of each 4 bits into those 4 bits 
    x = (x + (x >> 4)) & m4;        //put count of each 8 bits into those 8 bits 
    return (x * h01)>>56;  //returns left 8 bits of x + (x<<8) + (x<<16) + (x<<24) + ... 
}

popcount3进一步进行了优化，只看最后一步:return (x * h01)>>56;
此步之前的x=[a]8[b]8[c]8[d]8[e]8[f]8[g]8[h]8
x*h01 = x*0x0101010101010101 = x+(x<<8)+(x<<16)...+(x<<56)
即x=[a+b+c+d+e+f+g+h|L56], L56指低56位
右移56位获得a+b+c+d+e+f+h的值，得解。

此外还有比较有趣的算法：

//This is better when most bits in x are 0
//It uses 3 arithmetic operations and one comparison/branch per "1" bit in x.
int popcount_4(uint64 x) {
    int count;
    for (count=0; x; count++)
        x &= x-1;
    return count;
}

上面这个算法在已知0的数目比较多时候很高效。
此算法基于这样一个事实：x-1使得以二进制表示的x，从低向高位开始包括第一个1在内的值，都由0变成1，由1变成0。如11-01 = 10, 10 – 01 = 01, 01 – 01 = 00, 100 – 001 = 011。而&操作使得发生变化的位置都变成0，这样就去除了1个1，从而有几个1就&几次，最终x必变成0.
下面算法消除了popcount4的循环

//This is better if most bits in x are 0.
//It uses 2 arithmetic operations and one comparison/branch  per "1" bit in x.
//It is the same as the previous function, but with the loop unrolled.
#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;

最有趣的是查表法，当有足够的内存时，我们可以用空间换时间，从而得到O(1)的最优算法。
以4bit的串为例，可以构造一个数组int counts[16]={0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4}.
对于4bit的x,x的hamming weight即为：counts[x].
对于32bit的串，也可以使用分成两部分查表的方法来节省一点内存：

static unsigned char wordbits[65536] = { bitcounts of ints between 0 and 65535 };
static int popcount(uint32 i)
{
    return (wordbits[i&0xFFFF] + wordbits[i>>16]);
}

Hamming Weight还有很多应用，这里只是简单记录一下它在求解popcount上的用法。