Hamming Weight的算法分析(转载)

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

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

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

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


当然,还有更优的算法。

这个问题其实是HammingWeight的一个应用,又叫做populationcount,popcountorsidewayssumHammingWeight详见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语言实现。

[cpp] view plaincopy

    //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次算数操作,比起前面的算法来说要明显减少了。


[cpp] view plaincopy

    //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效果应该一样的吧..?)


[cpp] view plaincopy

    //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的值,得解。

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

[cpp] view plaincopy

    //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的循环


[cpp] view plaincopy

    //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上的用法。

posted @ 2015-03-18 19:43 禅意 阅读(...) 评论(...) 编辑 收藏