关于这个问题,网络上讨论的很多,可以找到大量的资料,我觉得就就是下面这一篇讲的最好,也非常的全面:
统计无符号整数二进制中 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