调整麻省理工学院的比特计数算法来并行计算单词数？

Question

我想使用众所周知的 MIT 位计数算法的一个版本，使用 SSE2 指令来计算 Conway 生命游戏中的邻居。

这是 c 中的 MIT 位计数，扩展为 count bitcounts > 63 位。

int bitCount(unsigned long long n)
{
unsigned long long uCount;

uCount = n – ((n >> 1) & 0×7777777777777777)
           - ((n >> 2) & 0×3333333333333333)
           - ((n >> 3) & 0×1111111111111111);
return ((uCount + (uCount >> 4))
& 0x0F0F0F0F0F0F0F0F) % 255;
}

这是 Pascal 的一个版本，

function bitcount(n: uint64): cardinal;
var ucount: uint64;
begin
  ucount:= n - ((n shr 1) and $7777777777777777)
             - ((n shr 2) and $3333333333333333) 
             - ((n shr 3) and $1111111111111111);
  Result:= ((ucount + (count shr 4)) 
           and $0F0F0F0F0F0F0F0F) mod 255;
end;

我希望并行计算该结构中的位。

  32-bit word where the pixels are laid out as follows.
  lo-byte         lo-byte neighbor
  0 4 8 C  048C   0 4 8 C 
   +---------------+
  1|5 9 D  159D   1|5 9 D 
   |               |
  2|6 A E  26AE   2|6 A E  
   +---------------+
  3 7 B F  37BF   3 7 B F 
 |-------------|            << slice A
   |---------------|        << slice B
     |---------------|      << slice C

请注意该结构中间有 16 位需要查找。 我想使用 SSE2 计算中间 16 位中每一位的邻居计数。 为了做到这一点，我将切片 A 放入 XMM0 低位双字中，将切片 B 放入 XXM0-dword1 等中。
我将 XMM0 复制到 XMM1，并屏蔽 XMM0 低位字中的位 5 的位 012-456-89A，对 XMM0 的 word1 执行相同操作，依此类推。不同的切片和掩码，以确保 XMM0 和 XMM1 中的每个单词保存不同像素的邻居。

问题
如何调整 MIT 位数以最终得到每个 XMM 字中每个字/像素的位数？

备注
我不想使用查找表，因为我已经有了这种方法，并且我想 测试 SSE2 是否会通过不需要对查找表进行内存访问来加速该过程。

使用 SSE 汇编的答案将是最佳的，因为我在 Delphi 中对此进行编程，因此我使用 x86+SSE2 汇编代码。

Answer 1

MIT 算法很难在 SSE2 中实现，因为没有可用于最终 ... % 255 表达式的整数模指令。在各种 popcnt 方法中，最适合 SSE 的方法可能是 Henry S. Warren 的“Hackers Delight”，我使用 SSE 内在函数在 C 中实现了它：

#include <stdio.h>
#include <emmintrin.h>

__m128i _mm_popcnt_epi16(__m128i v)
{
    v = _mm_add_epi16(_mm_and_si128(v, _mm_set1_epi16(0x5555)), _mm_and_si128(_mm_srli_epi16(v, 1), _mm_set1_epi16(0x5555)));
    v = _mm_add_epi16(_mm_and_si128(v, _mm_set1_epi16(0x3333)), _mm_and_si128(_mm_srli_epi16(v, 2), _mm_set1_epi16(0x3333)));
    v = _mm_add_epi16(_mm_and_si128(v, _mm_set1_epi16(0x0f0f)), _mm_and_si128(_mm_srli_epi16(v, 4), _mm_set1_epi16(0x0f0f)));
    v = _mm_add_epi16(_mm_and_si128(v, _mm_set1_epi16(0x00ff)), _mm_and_si128(_mm_srli_epi16(v, 8), _mm_set1_epi16(0x00ff)));
    return v;
}

int main(void)
{
    __m128i v0 = _mm_set_epi16(7, 6, 5, 4, 3, 2, 1, 0);
    __m128i v1;

    v1 = _mm_popcnt_epi16(v0);

    printf("v0 = %vhd\n", v0);
    printf("v1 = %vhd\n", v1);

    return 0;
}

编译和测试如下：

$ gcc -Wall -msse2 _mm_popcnt_epi16.c -o _mm_popcnt_epi16
$ ./_mm_popcnt_epi16 
v0 = 0 1 2 3 4 5 6 7
v1 = 0 1 1 2 1 2 2 3
$ 

看起来就像大约 16 条算术/逻辑指令一样，因此它应该以每点大约 16 / 8 = 2 个时钟运行。

如果需要，您可以轻松地将其转换为原始汇编程序 - 每个内在函数映射到单个指令。