popcount 算法分析

本文轉載自查看原文 2013-03-01 21:59 7619 algorithm

/* 轉載請注明出處： http://www.cnblogs.com/Martinium/articles/popcount.html */

簡介：　　　　

population count，簡稱 popcount 或叫 sideways sum. 是計算一個數的二進制表示有多少位是1，在一些場合下很有用，比如計算0-1稀疏矩陣(sparse matrix)或位數組(bit array)中非零元素個數、比如計算兩個字符串的漢明距離(Hamming distance)。然而 Intel 直到2008年11月 Nehalem 架構的處理器Core i7 才引入了 SSE4.2 指令集，其中就有 CRC32 和 POPCNT 指令，POPCNT可以處理16,32,64位整數。

關於popcount有很多算法，搜集到的大致有以下幾種：

iterated_popcnt
sparse_popcnt
dense_popcnt
lookup_popcnt
parallel_popcnt
nifty_popcnt
hacker_popcnt
hakmem_popcnt
assembly_popcnt

代碼分析：

int iterated_popcnt(unsigned int n)
{
    int count=0;
    for(;n;n>>=1)
        count+=n&1u;
    return count;
}

iterated_popcnt 通過迭代查找1累加，多少位就循環多少次，簡單明了也因此最慢。

int sparse_popcnt(unsigned int n)
{
    int count=0;
    while(n)
    {
        ++count;
        n&=(n-1);
    }
    return count;
}

sparse_popcnt 是對 iterated_popcnt 的改進，每次迭代總是將最右邊的非零位置零。減法的妙用。
試想一下，一個僅最高位為1的整數，用此方法的話僅需一次迭代；而 iterated_popcnt 還是會“乖乖的”迭代64次，對64位整形來講。

int dense_popcnt(unsigned int n)
{
    int count=CHAR_BIT*sizeof(unsigned int);
    n^=(unsigned int)-1;
    while(n)
    {
        --count;
        n&=(n-1);
    }
    return count;
}

dense_popcnt 只是 sparse_popcnt 的一個變種。
如果我們有先驗條件，知道數n的bit位有很多很多的1(取名 dense 的由來)，可以先取反求 popcount，最后用int位數一減就得到結果。如果數很隨機，即數的1、0位數差不多(not too sparse, not too dense)，將失去優勢，與 sparse_popcnt 表現一樣，如果忽略掉開始的一次取反操作。
n&(n-1) 式子還有一個用途，就是判斷一個整數 x 是否是 2 的 n 次冪。條件是 x!=0 && (x&(x-1))==0 。

int lookup_popcnt(unsigned int n)
{
#if 0 /* generate the table algorithmically, and you should put it outside. */
    static unsigned char TABLE[256]={0};
    for(unsigned int i=1;i<256;++i)
        TABLE[i]=TABLE[i>>1]+(i&1u);
    
    unsigned char *p=(unsigned char*)&n;
    return TABLE[p[0]]+TABLE[p[1]]+TABLE[p[2]]+TABLE[p[3]];
    
#else # define BIT2(n) n,            n+1,            n+1,            n+2 # define BIT4(n) BIT2(n), BIT2(n+1), BIT2(n+1), BIT2(n+2)
# define BIT6(n) BIT4(n), BIT4(n+1), BIT4(n+1), BIT4(n+2)
# define BIT8(n) BIT6(n), BIT6(n+1), BIT6(n+1), BIT6(n+2)

    static const unsigned char TABLE[256]={BIT8(0)};

    return TABLE[n       & 0xff]+
           TABLE[(n>>8)  & 0xff]+
           TABLE[(n>>16) & 0xff]+
           TABLE[(n>>24) & 0xff];
           
#endif          
}

lookup_popcnt 通過查表法，用空間換時間，所以來得快。
這里沒有制作整個 uint 區間的（需太大空間）表，而是一個 uchar 表，分四次求值然后相加。你可以嘗試 ushort 表，然后查表兩次(上面用了四次)，看會不會更快，畢竟直接索取結果應該比任何計算都要來得快，按理說。
上面的代碼列出了兩種方法構造 TABLE：
1. 利用 f(n)=f(n/2)+n%2, f(0)=0。直觀而且方便移植。不要求char一定為8位，即(CHAR_BIT==8)
2. 利用宏擴展遞歸展開，很漂亮的寫法。上面假定了char為8位。有些舊的機器一個字節7位，就需要作一些修改。
CHAR_BIT 和 UCHAR_MAX 作為宏定義在 limit.h 文件中，方便移植。

int lookup_popcnt(unsigned int n)
{
#define TBL_LEN (1<<CHAR_BIT)

#define BIT2(n)      n,       n+1,       n+1,       n+2
#define BIT4(n) BIT2(n), BIT2(n+1), BIT2(n+1), BIT2(n+2)
#define BIT6(n) BIT4(n), BIT4(n+1), BIT4(n+1), BIT4(n+2)
#define BIT8(n) BIT6(n), BIT6(n+1), BIT6(n+1), BIT6(n+2)

    static const unsigned char TABLE[TBL_LEN]=
        {
#if (CHAR_BIT==8)
            BIT8(0)
#elif (CHAR_BIT==7)
            BIT6(0), BIT6(1)
#else
# error "BITX to be added here."
#endif
        };

    return TABLE[n                 & UCHAR_MAX]+
           TABLE[(n>>CHAR_BIT)     & UCHAR_MAX]+
           TABLE[(n>>(CHAR_BIT*2)) & UCHAR_MAX]+
           TABLE[(n>>(CHAR_BIT*3)) & UCHAR_MAX];

#undef TBL_LEN
}

#define POW2(c) (1u<<(c))
#define MASK(c) (((unsigned int)(-1))/(POW2(POW2(c))+1u))
#define COUNT(x,c) ((x)&MASK(c)) + (((x)>>(POW2(c)))&MASK(c))

int parallel_popcnt(unsigned int n)
{
    n=COUNT(n,0);
    n=COUNT(n,1);
    n=COUNT(n,2);
    n=COUNT(n,3);
    n=COUNT(n,4);
/*  n=COUNT(n,5);  for 64-bit integers */
    return n ;
}

parallel_popcnt 利用並行計算。給出圖片展示：

想象一下已算好的每k位一組的 bitcount，向2*k位挺進。
num=a*2^k+b, a為高k位，b為低k位，於是結果為 a+b。
用關於num的表達式，用mask分別屏蔽高低位取出a和b再加。
如此，32位整數需5次迭代，64位整數需6次迭代。

關於第n次迭代MASK的值
[0]. 01010101 01010101 01010101 01010101 = 55555555 = FFFFFFFF/(0x2+1)
[1]. 00110011 00110011 00110011 00110011 = 33333333 = FFFFFFFF/(0x4+1)
[2]. 00001111 00001111 00001111 00001111 = 0F0F0F0F = FFFFFFFF/(0x10+1)
[3]. 00000000 11111111 00000000 11111111 = 00FF00FF = FFFFFFFF/(0x100+1)
[4]. 00000000 00000000 11111111 11111111 = 0000FFFF = FFFFFFFF/(0x10000+1)
上面除了左邊第一列是二進制外，右邊的都是十六進制數。
第n次的數的規律是 2ⁿ個0，接着 2ⁿ個1，然后循環。2^2^n對應高2ⁿ位，1對應低2ⁿ位。
所以右邊的除數是 2^2^n+1，此處^為乘方，右結合。其實就是費馬數(Fermat number)。

#define MASK_01010101 (((unsigned int)(-1))/3)
#define MASK_00110011 (((unsigned int)(-1))/5)
#define MASK_00001111 (((unsigned int)(-1))/17)

int nifty_popcnt(unsigned int n)
{
    n = (n & MASK_01010101) + ((n>>1) & MASK_01010101) ;
    n = (n & MASK_00110011) + ((n>>2) & MASK_00110011) ;
    n = (n & MASK_00001111) + ((n>>4) & MASK_00001111) ;        
    return n%255 ;
}

nifty_popcnt 與 parallel_popcnt 開始都一樣，只是寫法不同，用3次迭代完成了每8位一組的和，不同的在於末尾用模255一下結束了運算，較之漂亮(nifty)。這里面涉及到一些數學知識。
K進制數B (B_nB_n-1...B₁B₀)，B=B_n*K^n-1+B_n-1*K^n-2+...+B₁*K+B₀Kⁱ≡1(mod K-1),可以利用二項式Kⁱ=((K-1)+1)ⁱ展開，或用數學歸納法證明此結論。
B (mod K-1)
=B_n*K^n-1+B_n-1*K^n-2+...+B₁*K+B₀(mod K-1)
=B_n+B_n-1+...+B₁+B₀(mod K-1)
結論：一個K進制數模(K-1)的結果，等於此K進制數的各位相加再模K-1
於是，popcnt(B)≡B (mod K-1)
當然，我們能肯定 B_n+B_n-1+...+B₁+B₀<K-1 的話，可以加強結論：popcnt(B)=B%(K-1)
上面每8位一組，相當於2⁸=256進制，所以用了255這個數；為了使用上面的等式計算，必須至少得3次迭代。
2次迭代創造 2⁴=16 進制，而對於一個32位整形，popcount 最大值為32，32>16-1；32<256-1
想想一下64位整形，popcount 可能的最大取值是64，這里要選取的數是511。

int hacker_popcnt(unsigned int n)
{
    n -= (n>>1) & 0x55555555;
    n  = (n & 0x33333333) + ((n>>2) & 0x33333333);
    n  = ((n>>4) + n) & 0x0F0F0F0F;
    n += n>>8;
    n += n>>16;
    return n&0x0000003F;
}

取名 hacker_popcnt 是因為來自於 Hacker's Delight 這本書。網上有電子版下載，強烈推薦翻閱。
書中第5章專門講 Counting Bits。上面的那段代碼來自第一節的 Figure 5-2 Counting 1-bits in a word.
該代碼也被 Android 開源項目所用，路徑 <AOSP>libcore/luni/src/main/java/java/lang/
FFmpeg 中也有這段代碼，路徑 libavutil/common.h 可以 grep 一下關鍵字 popcount。
其中的 Integer 和 Long 類的方法 bitCount，分別對應32位和64位的 popcount
第一行的減法的由來：二進制二位數 2a+b 左移一位得 a，於是 (2a+b)-a=a+b
方法原理同上面的 nifty_popcnt 一樣，這里不多描述。

/* HAKMEM Popcount

  Consider a 3 bit number as being
        4a+2b+c
  if we shift it right 1 bit, we have
        2a+b
  subtracting this from the original gives
        2a+b+c
  if we shift the original 2 bits right we get
        a
  and so with another subtraction we have
        a+b+c
  which is the number of bits in the original number.

  Suitable masking  allows the sums of  the octal digits  in a 32 bit  number to
  appear in  each octal digit.  This  isn't much help  unless we can get  all of
  them summed together.   This can be done by modulo  arithmetic (sum the digits
  in a number by  molulo the base of the number minus  one) the old "casting out
  nines" trick  they taught  in school before  calculators were  invented.  Now,
  using mod 7 wont help us, because our number will very likely have more than 7
  bits set.   So add  the octal digits  together to  get base64 digits,  and use
  modulo 63.   (Those of you  with 64  bit machines need  to add 3  octal digits
  together to get base512 digits, and use mod 511.)
 
  This is HACKMEM 169, as used in X11 sources.
  Source: MIT AI Lab memo, late 1970's.
*/

int hakmem_popcount(unsigned int n)
{
    unsigned int tmp;
    tmp = n - ((n>>1) & 033333333333) - ((n>>2) & 011111111111);
    return ((tmp + (tmp>>3)) & 030707070707) % 63;
}

取名 hakmem_popcnt 是因為來自於 HAKMEM (Hacks Memo)，來自 MIT AI Lab，里面有很多關於數學的有用而且漂亮的算法。你可以多回味一下這個算法。^_^
這里用到了一個技巧：popcnt(n) = n - (n>>1) - (n>>2) - ... - (n>>31) = a₃₁+a₃₀+...+a₀
證明不難，提示：2^k-2^k-1-2^k-2-...-2¹-2⁰=1，對任意整數k成立。
上面先 triple 一下 popcnt(abc)=popcnt(4a+2b+c)=(4a+2b+c)-(2a+b)-a=a+b+c，接着再 double 一下，最后模63。
從上面的 nifty_popcnt 分析得知，其實 popcount(unsigned int)<=32，這里假設 unsigned int是32位整形。用63這個
數模就夠了，因為63是滿足2ⁿ-1的最小數。於是n=6，表示我們需要獲得每6位一求的結果。6=3x2，這就是上面先 triple
再 double 的原因。而 nifty_popcnt 總是 double (x2)。

int assembly_popcnt(unsigned int n)
{
#ifdef _MSC_VER /* Intel style assembly */
    __asm popcnt eax,n;
#elif __GNUC__ /* AT&T style assembly */
    __asm__("popcnt %0,%%eax"::"r"(n));
}

最終來到了匯編指令 POPCNT，匯編語言有兩種style：Intel 和 AT&T。語法一樣，只是寫法不同。之間的區別自己 Google。
Visual C++ 用的是 Intel style；類UNIX 使用 AT&T style。大學學的是前面一種。個人覺得 Intel 的簡單易懂，上面的一行指令足已證明。C/C++函數的返回值都是存入 EAX 寄存器的。函數很簡短，所以被調用的時候可以直接 inline 到其他函數。

關於匯編指令的用法，請查閱 Intel 官方手冊，地址在這里。你的機器不一定支持（編譯通過但運行出錯），所以最好先檢測一下。 Linux 平台通過 cat /proc/cpuinfo | grep popcnt 命令就可以看到是否支持。Before an application attempts to use the POPCNT instruction, it must check that the processor supports SSE4.2 (if CPUID.01H:ECX.SSE4_2[bit 20] = 1) and POPCNT (if CPUID.01H:ECX.POPCNT[bit 23] = 1).

asm 本來是 C++ 的關鍵字的，MSDN上說 The __asm keyword replaces C++ asm syntax. asm is reserved for compatibility with other C++ implementations, but not implemented. Use__asm. 如果用 GCC 編譯器的話，可以使用asm，但編譯的時候要用 gcc popcnt.c -o bitcnt -std=c99 -fasm，-fasm是讓編譯器認"asm", "inline" or "typeof"為關鍵字。
GNU編譯器也內置了很多函數，也包括 int __builtin_popcount (unsigned int x);，如果你的機器架構支持的話，會直接譯成一條CPU指令。

測試:
上面的代碼放在一起做個測試，比較一下效率。
僅測試32位整形(unsigned int)，你可以修改測試一下64位整形(unsigned long long int)。Microsoft Visual C++的64位編譯器不支持內嵌匯編，可以通過編譯器提供的 intrinsic function 或者用32位分兩次求和。
注意先驗證方法的正確性再做測試。使用 CLOCKS_PER_SEC 宏，因為

Linux 平台 /usr/include/bits/time.h #define CLOCKS_PER_SEC 1000000l
Windows 平台 Microsoft Visual Studio 10.0\VC\include\time.h #define CLOCKS_PER_SEC 1000

完整的程序從這里下載，或者在線閱讀 popcount.cpp

編譯運行的結果：

總結：
對比 popcount 的各種算法，高效在於能利用並行計算，去掉循環，使用減法和模運算。
通過減1的循環算法(parse/dense)在知道數只有三五位為1（或0）的話，其實效率也不賴。
查表法的效率挺不錯的，如果沒有硬件指令的支持，選用這個是可以的。
Hacker's Delight 中的算法，在開源項目中廣為引用。
算法就要充分挖掘里面的數學知識。

參考：

The Aggregate Magic Algorithms by enry Gordon Dietz at University of Kentucky
Puzzle: Fast Bit Counting Delightful Puzzles
HAKMEM (bit counting is memo number 169), MIT AI Lab, Artificial Intelligence Memo No. 239, February 29, 1972.
Bit Twiddling Hacks by Sean Anderson at Stanford University.
Bitwise Tricks and Techniques by Donald E. Knuth (The Art of Computer Programming, Part IV).
Hacker's Delight Hacker's_Delight_2nd.pdf

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 算法與算法分析 join算法分析算法分析 offsetParent算法分析 TimSort算法分析眾數的算法分析共識算法分析算法分析-漸進分析自動排課算法分析 Cuckoo hash算法分析