Sunday, November 29, 2009

Fast bit counting

最近 mesa3d-dev mailing list 有一個討論串在討論 bit counting,也就是在算一個 unsigned int 裡有幾個 bit 是 1。

假設 unsigned int 有 32 bits 寬,最簡單的算法就是一個迴圈跑 32 次,一個一個 bit 去看。想當然這不是有效率的算法。事實上,只要做簡單的修改,就可以讓迴圈跑的次數降到跟 bit 為 1 的 bit 數相同。換句話說,如果 unsigned int 裡只有 3 個 bit 為 1,迴圈只要跑三次就可以。不過像這麼基礎的運算,一定有不少人下過功夫在找最快的算法。一個問題,工程師會去尋找更好的方法;但是身為鄉民,我們感興趣的只有最好的方法,而且我們從不自己找!很幸運地,在最開頭提到的討論串裡頭,就有人提出一個號稱最快的算法。

討論串中提到的方法經過簡化可以寫成

int bitcount(unsigned int v)
{
v = v - ((v >> 1) & 0x55555555);
v = (v & 0x33333333) + ((v >> 2) & 0x33333333);
v = (v + (v >> 4)) & 0x0f0f0f0f;
return (v * 0x01010101) >> 24;
}


在深入去看這段程式碼前,最好可以用抽像點的方式看它想做的事



上圖中,v0 是使用者的輸入。展開成二進位表示法,每個 bit 的值不是 1 就是 0。如果有一個方法可以讓所有的 bit 變成兩個兩個一組 (v1)、4 個 4 個一組 (v2)、8 個 8 個一組 (v3)、一直到 32 個 bit 自己成一組 (v5),那 v5 的值剛好就是我們要的結果。從 v4 看起,假設我們已經有 v4 了,那麼應該不難發現
v5 = (v4 & 0xffff) + (v4 >> 16);


以此類推,很快就可以發現我們需要的是

int bitcount2(unsigned int v)
{
v = (v & 0x55555555) + ((v >> 1) & 0x55555555); /* v1 */
v = (v & 0x33333333) + ((v >> 2) & 0x33333333); /* v2 */
v = (v & 0x0f0f0f0f) + ((v >> 4) & 0x0f0f0f0f); /* v3 */
v = (v & 0x00ff00ff) + ((v >> 8) & 0x00ff00ff); /* v4 */
v = (v & 0x0000ffff) + ((v >> 16) & 0x0000ffff); /* v5 */
return v;
}


這時再回過頭去看 bitcount,可以注意到雖然算法不完全相同,但它前三行在做的也是求 v1、v2 與 v3。但它接下來不是求 v4,而是把 v3 乘上 0x01010101。透過四則運算,這個乘法可以改寫成加法

v = (v3 << 24) + (v3 << 16) + (v3 << 8) + v3


把 v3 的展開式代入,會發現



將這個結果向右 shift 24 bits 得到的也是 v5。這就是 bitcount 的做法。

如果把寫死的 magic number 換成適當型別的除法與補數運算,這個 bitcount 函數對 32/64/128-bits 寬的 unsigned int 都可以適用。寬度的限制是來自於 v4' 的第一個小括號,它不能大於或等於 2^8。不過如果 bitcount 求到 v4 再來做乘法,這個限制也可以被放寬。