有效地交错位答案

【问题标题】：Interleave bits efficiently有效地交错位
【发布时间】：2022-05-17 12:17:31
【问题描述】：

我需要在 2 个uint32_t 中生成uint64_t 交错位：如果A=a0a1a2...a31 和B=b0b1...b31，我需要C=a0b0a1b1...a31b31。有没有办法有效地做到这一点？到目前为止，我只有一个幼稚的方法，其中包含 32 次迭代的 for 循环，其中每次迭代都执行 C|=((A&(1<<i))<<i)|((B&(1<<i))<<(i+1))。

我想应该有一些数学技巧，例如将 A 和 B 乘以某个特殊数字，这会导致在生成的 64 位数字中将它们的位与零交错，因此剩下的就是 or 这些产品。但是我找不到这样的乘数。

另一种可能的方法是编译器内部指令或汇编指令，但我不知道。

【问题讨论】：

This 看起来可能有效。不确定它是否或如何扩展。
它应该很容易扩展；添加 16 的初始移位和 0x0000ffff0000ffff 的掩码，并适当地扩展所有其他掩码。因此，新的第一步将高 16 位分割并移动到最高字的低位。
可以使用x86的BMI2指令集吗？两个 64 位 PDEP with opposite masks 和一个普通的按位或应该可以解决问题。 PDEP 的内在函数是 _pdep_u64()
@PeterCordes，这是一个出色的理论解决方案，但实际上我在 Phenom 2 处理器和 CUDA 中都没有 BMI2/PDEP。

标签： c++ algorithm math assembly bit-manipulation

【解决方案1】：

NathanOliver 的链接提供了 16 位 -> 32 位实现：

static const unsigned int B[] = {0x55555555, 0x33333333, 0x0F0F0F0F, 0x00FF00FF};
static const unsigned int S[] = {1, 2, 4, 8};

unsigned int x; // Interleave lower 16 bits of x and y, so the bits of x
unsigned int y; // are in the even positions and bits from y in the odd;
unsigned int z; // z gets the resulting 32-bit Morton Number.  
                // x and y must initially be less than 65536.

x = (x | (x << S[3])) & B[3];
x = (x | (x << S[2])) & B[2];
x = (x | (x << S[1])) & B[1];
x = (x | (x << S[0])) & B[0];

y = [the same thing on y]

z = x | (y << 1);

其工作原理：

将 x 的低 8 位保留在原处。将高 8 位上移 8 位；
分成两半，做同样的事情，这次将 4 位的低位对留在原处，将其他位向上移动 4 位；
一次又一次。

即它进行如下：

   0000 0000 0000 0000  abcd efgh ijkl mnop
-> 0000 0000 abcd efgh  0000 0000 ijkl mnop
-> 0000 abcd 0000 efgh  0000 ijkl 0000 mnop
-> 00ab 00cd 00ef 00gh  00ij 00kl 00mn 00op
-> 0a0b 0c0d 0e0f 0g0h  0i0j 0k0l 0m0n 0o0p

然后将两个输入组合在一起。

根据我之前的评论，要将其扩展到 64 位，只需添加 16 的初始移位和 0x0000ffff0000ffff 的掩码，因为您可以直观地遵循模式或作为分而治之的步骤，转动将 32 位问题分成两个不重叠的 16 位问题，然后使用 16 位解决方案。

【讨论】：

分而治之方法使所需操作的数量翻倍，而使用第 5 个数组项和or-and-assign 进行扩展，如果有效，则将操作数量增加 1.25 倍。
在这种情况下，分而治之不会翻倍，因为 32 位值被分成两个 16 位量，适合整个 64 位字，您的 CPU 可以原子地对其进行操作。但我们在争论语义。我把它算作逻辑上的一个初始除法，即使被除法的东西仍然适合一个标量。
谢谢，我一开始误解了。我虽然您的意思是先在低 32 位上执行这 4 个 or-and-assign 操作，然后在高 32 位上执行。
@J.Schultke 该问题已被标记为[c++]，因此这是未标记块的默认代码格式，没有 cmets 或反引号语言。我认为您的编辑根本没有改变第一个代码块的语法突出显示。其他小的更改是改进，因此不会回滚，但将来您可以通过牢记这一点并寻找实际的不良突出显示来减少您的编辑噪音。

【解决方案2】：

对于较大的整数，值得一提的是用于有限域乘法（无进位乘法）的 clmul x86 扩展。将整数与零交错等效于整数与自身的无进位乘法，这是一条 ALU 指令。

【讨论】：

不错的选择，64 位输入/128 位输出甚至非常理想。但同样值得注意的是，这目前仍然是一个有点“新”的扩展，并且不会在您可能仍然在野外遇到的 AMD Phenom II 和 Intels Ivy Bridge 处理器上提供。
@Ext3h: felixcloutier.com/x86/pclmulqdq 自 Westmere for Intel（第二代 Nehalem）起可用。 godbolt.org/z/jWaGWv。但是，是的，它是一个扩展，而不是 x86-64 基线。
由于我们提到了 x86 扩展，BMI2 pdep r64, r64, r/m64 在 Intel（1 uop，3c 延迟）上有效，但在 AMD 上无效。它可以将 32 位扩展为 64 位整数，其间有零。（BMI2 在 Haswell 及更高版本以及 AMD Zen 及更高版本上可用，但 pdep/pext 在 Zen 上非常慢）

【解决方案3】：

一个简短的、预先计算好的数组查找算作“数学技巧”吗？

预先计算一个包含 256 个uint16_ts 的数组：

static const uint16_t lookup[256]={0x0000, 0x0001, 0x0005 ..., 0x5555};

我们可以将两个 8 位的值交错，很容易得出一个 16 位的值：

uint16_t interleave(uint8_t a, uint8_t b)
{
    return (lookup[a] << 1) | lookup[b];
}

如何扩展它以将两个 32 位值交织成一个 64 位值应该很明显：调用此方法四次，对于构成 uint32_t 的四个字节中的每一个，然后是 << 和 @ 987654326@结果一起。贿赂编译器以内联整个事情，最终结果应该是相当快速和便宜的。

由于现在 RAM 很便宜，您可能还需要考虑预先计算的 65536 uint32_ts 表。

【讨论】：

内存很便宜；缓存未命中不是。大多数程序区域的内存瓶颈已经超过了 CPU。这可能会赢得一个微基准测试，其中没有其他代码会占用您的缓存空间，但它不会在现实生活中获胜，除非此操作在一个非常紧密的循环中以突发方式发生。
在某些情况下仍然是个好主意。我需要交错两个字节，但通过对每个字节的上下 4 位半字节进行 2 次查找可能是值得的。
@MadsY：确实，4 位 LUT（16 个 1 字节条目）足够小，可以在缓存中可靠地保持热状态，并且可以在内存和计算之间进行很好的权衡（以及 + 转移到隔离每个半字节，或者合并）。

【解决方案4】：

为了使 saolof 的回答具体化，以下是使用 CLMUL 指令集的实现，每次调用交错两对 uint32_ts：

#include <immintrin.h>
#include <stdint.h>

typedef struct {
  uint32_t x;
  uint32_t y;
} uint32_2;

static inline void interleave_clmul(uint32_2 *input, uint64_t *out) {
  __m128i xy = _mm_load_si128((const __m128i *)input);

  xy = _mm_shuffle_epi32(xy, 0b11011000);

  // Two carryless multiplies
  __m128i p2 = _mm_clmulepi64_si128(xy, xy, 0x11);
  __m128i p1 = _mm_clmulepi64_si128(xy, xy, 0x00);

  // Bitwise interleave the results
  p2 = _mm_slli_epi16(p2, 1);
  __m128i p3 = _mm_or_si128(p1, p2);

  _mm_storeu_si128((__m128i *)out, p3);
}

编译为the following:

interleave_clmul(uint32_2*, unsigned long*):
        vpshufd         xmm0, xmmword ptr [rdi], 216    # xmm0 = mem[0,2,1,3]
        vpclmulqdq      xmm1, xmm0, xmm0, 17
        vpclmulqdq      xmm0, xmm0, xmm0, 0
        vpaddw          xmm1, xmm1, xmm1
        vpor            xmm0, xmm0, xmm1
        vmovdqu         xmmword ptr [rsi], xmm0
        ret

正如简单性所暗示的那样，在我的（Haswell）机器上，这个系统比使用pdep 指令的相应实现要快得多。时序非常嘈杂，并且对所涉及的阵列的大小很敏感，但是对于 10000 对的阵列，它是一致的 2.4 倍加速（pdep 为每对 2.27 个周期，此 impl 为 0.9 个周期）。

我使用了在 Daniel Lemire's blog 上找到的测试工具。

【讨论】：