[快速閱讀六] 統計記憶體資料中二進位制1的個數（SSE指令集最佳化版）.

　　關於這個問題，網路上討論的很多，可以找到大量的資料，我覺得就就是下面這一篇講的最好,也非常的全面：

　 統計無符號整數二進位制中 1 的個數（Hamming Weight）

　　在指令集不參與的情況下，分治法速度最快，MIT HAKMEM 169 演算法因為最後有一個mod取餘操作，速度要稍微慢一點，256元素的查表演算法速度要次之，當然，其實要建議那個256元素的表不要使用int型別，而是使用unsigned char型別的，儘量減少表的記憶體佔用量，也意味著cache miss小一些。 16位的查表演算法速度反而慢了不少，主要是因為他用while，即使我們把他展開，也需要8次資料組合，還是比16位的慢。其他的就不要說了，都比較慢。

　　在SSE4指令集能得到CPU的支援時，可以有一個直接的指令_mm_popcnt_u32可以使用，這個就可以加速很多了，一個常用的過程如下：

    Amount = 0;
    for (int Y = 0; Y < Length; Y++)
    {
        Amount += _mm_popcnt_u32(Array[Y]);
    }

　　一千萬的隨機資料，用這個指令大概只要3.2毫秒多就可以處理完成，如果稍微改下程式碼，讓他能並行化一點，如下所示：

    Amount = 0;
    for (int Y = 0; Y < Length; Y+=4)
    {
        int T0 = _mm_popcnt_u32(Array[Y + 0]);
        int T1 = _mm_popcnt_u32(Array[Y + 1]);
        int T2 = _mm_popcnt_u32(Array[Y + 2]);
        int T3 = _mm_popcnt_u32(Array[Y + 3]);
        Amount += T0 + T1 + T2 + T3;
    }

　　還可以提高到大約2.7ms就可以處理完。

　　其實，現在在執行的新的CPU基本上沒有那個不支援SSE4的了，但是也不排除還有一些老爺機。因為SSE4最早是2008年釋出的，如果CPU不支援SSE4，但是支援SSE3(2004年釋出的），那是否有合適的指令集能加速這個過程呢，實際上也是有的。

　　我們這裡喵上了統計無符號整數二進位制中 1 的個數（Hamming Weight）一文中的16元素查表演算法，原文中的程式碼為：

　　Amount = 0;
    for (int Y = 0; Y < Length; Y++)
    {
        unsigned int Value = Array[Y];
        while (Value)
        {
            Amount += Table16[Value & 0xf];
            Value >>= 4;
        }
    }

　　這個明顯是不合適指令處理的，前面說了，這個可以展開，展開後形式如下：

　　Amount = 0;
    for (int Y = 0; Y < Length; Y++)
    {
        unsigned int Value = Array[Y];
        Amount += Table16[Value & 0xf] + Table16[(Value >> 4) & 0xf] + Table16[(Value >> 8) & 0xf] + Table16[(Value >> 12) & 0xf] +
                Table16[(Value >> 16) & 0xf] + Table16[(Value >> 20) & 0xf] + Table16[(Value >> 24) & 0xf] + Table16[(Value >> 28) & 0xf];
    }

　　仔細觀察他的意思就是提取記憶體的4位，然後根據4位的值來查16個元素的表，我在之前的多個文章裡都有提高，16個元素的表（表內的值不能超過255）是可以借用一個_mm_shuffle_epi8指令進行最佳化的，一次性得到16個值。

　　_mm_shuffle_epi8是在SSE3就開始支援的，因此，這一改動可以將硬體的適應性提前4年。

　　具體的來首，就是我們載入16個位元組資料，然後和0xF進行and操作，得到每個位元組的低4位，然後進行shuffle，得到每個位元組低4位的二進位制中1的個數，然後在把原始位元組數右移4位，再和0xF進行and操作，得到每個位元組的高4位，然後進行shuffle，兩次shuffle的結果相加，就得到了這16個位元組資料的二進位制中1的個數。 具體程式碼如下所示：

    __m128i Table = _mm_loadu_si128((__m128i*)Table16);
    __m128i Mask = _mm_set1_epi8(0xf);
    __m128i UsedV = _mm_setzero_si128();
    for (int Y = 0; Y < Length; Y += 4)
    {
        __m128i Src = _mm_loadu_si128((__m128i*)(Array + Y));
        __m128i SrcL = _mm_and_si128(Src, Mask);
        __m128i SrcH = _mm_and_si128(_mm_srli_epi16(Src, 4), Mask);
        __m128i ValidL = _mm_shuffle_epi8(Table, SrcL);
        __m128i ValidH = _mm_shuffle_epi8(Table, SrcH);
        UsedV = _mm_add_epi32(UsedV, _mm_add_epi32(_mm_sad_epu8(ValidL, _mm_setzero_si128()), _mm_sad_epu8(ValidH, _mm_setzero_si128())));
    }
    //    提取出前面的使用SSE指令計算出的總的有效點數
    Amount = _mm_cvtsi128_si32(_mm_add_epi32(UsedV, _mm_unpackhi_epi64(UsedV, UsedV)));

　　注意到這裡的函式除了_mm_shuffle_epi8，其他的都是SSE2就已經能支援的了，其中_mm_sad_epu8可以快速的把16位元組結果相加。

　　 使用這個程式碼，測試上述1千萬資料，大概只需要2.1ms就能處理完，比最佳化後的_mm_popcnt_u32還要快。

　　 實際上，我還遇到一種情況，一個AMD的早期CPU，用CPUID看他支援的指令集，他是支援SSE4.2的，也支援SSE3，但是執行_mm_shuffle_epi8確提示不識別的指令，也是很奇怪，或者說如果遇到一個機器不支援SSE3,只支援SSE2，那是否還能用指令集最佳化這個演算法呢（SSE2是2001年釋出的）。

　　其實也是可以的，我們觀察上面的不使用指令集的版本的，特別是分冶法的程式碼：

　　Amount = 0;
    for (int Y = 0; Y < Length; Y++)
    {
        unsigned int Value = Array[Y];
        Value = (Value & 0x55555555) + ((Value >> 1) & 0x55555555);
        Value = (Value & 0x33333333) + ((Value >> 2) & 0x33333333);
        Value = (Value & 0x0f0f0f0f) + ((Value >> 4) & 0x0f0f0f0f);
        Value = (Value & 0x00ff00ff) + ((Value >> 8) & 0x00ff00ff);
        Value = (Value & 0x0000ffff) + ((Value >> 16) & 0x0000ffff);
        Amount += Value;
    }

　　這裡就是簡單的一些位運算和移位，他們對應的指令集在SSE2裡都能得到支援的，而且這個改為指令集也是水到渠成的事情：

　　UsedV = _mm_setzero_si128();
    for (int Y = 0; Y < Length; Y += 4)
    {
        __m128i Value = _mm_loadu_si128((__m128i*)(Array + Y));
        Value = _mm_add_epi32(_mm_and_si128(Value, _mm_set1_epi32(0x55555555)), _mm_and_si128(_mm_srli_epi32(Value, 1), _mm_set1_epi32(0x55555555)));
        Value = _mm_add_epi32(_mm_and_si128(Value, _mm_set1_epi32(0x33333333)), _mm_and_si128(_mm_srli_epi32(Value, 2), _mm_set1_epi32(0x33333333)));
        Value = _mm_add_epi32(_mm_and_si128(Value, _mm_set1_epi32(0x0f0f0f0f)), _mm_and_si128(_mm_srli_epi32(Value, 4), _mm_set1_epi32(0x0f0f0f0f)));
        Value = _mm_add_epi32(_mm_and_si128(Value, _mm_set1_epi32(0x00ff00ff)), _mm_and_si128(_mm_srli_epi32(Value, 8), _mm_set1_epi32(0x00ff00ff)));
        Value = _mm_add_epi32(_mm_and_si128(Value, _mm_set1_epi32(0x0000ffff)), _mm_and_si128(_mm_srli_epi32(Value, 16), _mm_set1_epi32(0x0000ffff)));
        UsedV = _mm_add_epi32(UsedV, Value);
    }
    //    提取出前面的使用SSE指令計算出的總的有效點數
    //Amount = _mm_cvtsi128_si32(_mm_add_epi32(UsedV, _mm_unpackhi_epi64(UsedV, UsedV)));
    Amount = UsedV.m128i_u32[0] + UsedV.m128i_u32[1] + UsedV.m128i_u32[2] + UsedV.m128i_u32[3];

　　這裡唯一要注意的就是最後從UsedV 變數得到Amount的過程不能用之前shuffle那一套程式碼了，因為這裡的UsedV的最高位裡含有了符號位，所以要換成下面哪一種搞法。

　　我們注意到，編譯執行這個程式碼後，我們得到的耗時大概是5.2ms，但是同樣的資料，前面的分冶法對應的C程式碼也差不多是5.5ms左右，速度感覺毫無提高，這是怎麼回事呢，我們嘗試反彙編C的程式碼，結果發現如下片段：

　　　　　　　　

　　注意到pand psrld paddd等指令沒有，那些就對應了_mm_and_si128 _mm_srli_epi32 _mm_add_epi32等等，原來是編譯器已經幫我們向量化了，而且即使在設定裡設定 無增強指令 (/arch:IA32) 選項，編譯器也會進行向量化（VS2019）。

　　

　　所以我暫時還得不到這個和純C比的真正的加速比。

　　但是，在編譯器沒有這個向量化能力時，直接手工嵌入SSE2的指令，還是能有明顯的加速作用的，不過也可以看到，SSE2的最佳化速度還是比SSE3的shuffle版本慢一倍的，而sse3的shuffle確可以比SSE4的popcount快30%。

　　以前我一直在想，這個演算法有什麼實際的應用呢，有什麼地方我會用到統計二進位制中1的個數呢，最近確實遇到過了一次。

　　具體的應用是，我有一堆資料，我要統計出資料中符合某個條件（有可能是多個條件）的目標有多少個，這個時候我們多次應用了_mm_cmpxx_ps等函式組合，最後得到一個Mask，這個時候我們使用_mm_movemask_ps來得到一個標記，我們看看_mm_movemask_ps 這個函式的具體意思：

　　他返回的是一個0到15之間的整形資料，很明顯我們可以把他儲存到一個unsigned char型別的變數裡，這樣，在計算完一堆資料後，我們就得到了一個mask陣列，這個時候我們統計下陣列裡有多少個二進位制1就可以得到符合條件的目標數量了。 當然，這裡和前面的還不太一樣，這個usnigned char變數的高4位一直為0，還可以不用處理的，還能進一步加速。

　　當然，如果系統支援AVX2，那還可以進一步做速度最佳化。這個就不多言了。

　　最後，列一下各個演算法的耗時比較資料吧：

　　相關測試程式碼地址: 資料流二進位制中1的個數統計

如果想時刻關注本人的最新文章，也可關注公眾號或者新增本人微信： laviewpbt

翻譯

搜尋

複製