比较两对 4 个变量并返回匹配数？

Question

给定以下结构：

struct four_points {
    uint32_t a, b, c, d;
}

比较两个这样的结构并返回匹配的变量数量（在任何位置）的绝对最快方法是什么？

例如：

four_points s1 = {0, 1, 2, 3};
four_points s2 = {1, 2, 3, 4};

我会寻找 3 的结果，因为三个数字在两个结构之间匹配。但是，鉴于以下情况：

four_points s1 = {1, 0, 2, 0};
four_points s2 = {0, 1, 9, 7};

那么我希望结果只有 2，因为在任一结构之间只有两个变量匹配（尽管第一个有两个零）。

我已经想出了一些用于执行比较的基本系统，但这是在短时间内被称为数百万次并且需要相对较快的系统。我目前最好的尝试是使用排序网络对任一输入的所有四个值进行排序，然后遍历排序后的值并保持相等的值的计数，从而相应地推进任一输入的当前索引。

是否有任何一种技术可能比排序和迭代性能更好？

Answer 1

在现代 CPU 上，有时正确应用蛮力是可行的方法。诀窍是编写不受指令延迟限制的代码，仅受吞吐量限制。

---

重复是否常见？如果它们非常罕见，或者有模式，使用分支来处理它们会使常见情况更快。如果它们真的不可预测，最好做一些无分支的事情。我正在考虑使用分支来检查它们罕见的位置之间的重复，并在更常见的位置使用无分支。

基准测试很棘手，因为带有分支的版本在用相同的数据测试一百万次时会大放异彩，但在实际使用中会有很多分支错误预测。

---

我还没有对任何东西进行基准测试，但是我想出了一个版本，它通过使用 OR 而不是加法来跳过重复项 来组合找到的匹配项。它编译为 gcc 完全展开的漂亮 x86 asm。 （没有条件分支，甚至没有循环）。

Here it is on godbolt。 （g++ 是愚蠢的，并且在 x86 setcc 的输出上使用 32 位操作，它只设置低 8 位。这种部分寄存器访问会降低速度。而且我什至不确定它是否会将高 24 位归零。 .. 无论如何，来自 gcc 4.9.2 的代码看起来不错，godbolt 上也是如此）

// 8-bit types used because x86's setcc instruction only sets the low 8 of a register
// leaving the other bits unmodified.
// Doing a 32bit add from that creates a partial register slowdown on Intel P6 and Sandybridge CPU families
// Also, compilers like to insert movzx (zero-extend) instructions
// because I guess they don't realize the previous high bits are all zero.
// (Or they're tuning for pre-sandybridge Intel, where the stall is worse than SnB inserting the extra uop itself).

// The return type is 8bit because otherwise clang decides it should generate
// things as 32bit in the first place, and does zero-extension -> 32bit adds.
int8_t match4_ordups(const four_points *s1struct, const four_points *s2struct)
{
    const int32_t *s1 = &s1struct->a; // TODO: check if this breaks aliasing rules
    const int32_t *s2 = &s2struct->a;
    // ignore duplicates by combining with OR instead of addition
    int8_t matches = 0;

    for (int j=0 ; j<4 ; j++) {
        matches |= (s1[0] == s2[j]);
    }

    for (int i=1; i<4; i++) { // i=0 iteration is broken out above
        uint32_t s1i = s1[i];

        int8_t notdup = 1; // is s1[i] a duplicate of s1[0.. i-1]?
        for (int j=0 ; j<i ; j++) {
            notdup &= (uint8_t) (s1i != s1[j]);  // like dup |= (s1i == s1[j]); but saves a NOT
        }

        int8_t mi = // match this iteration?
            (s1i == s2[0]) |
            (s1i == s2[1]) |
            (s1i == s2[2]) |
            (s1i == s2[3]);
    // gcc and clang insist on doing 3 dependent OR insns regardless of parens, not that it matters

        matches += mi & notdup;
    }
    return matches;
}

// see the godbolt link for a main() simple test harness.

---

在具有 128b 个向量的机器上，可以处理 4 个压缩的 32 位整数（例如 x86 和 SSE2），您可以将 s1 的每个元素广播到它自己的向量，去重，然后进行 4 次压缩比较。 icc 做了这样的事情来自动矢量化我的 match4_ordups 函数（在 Godbolt 上查看。）

使用 movemask 将比较结果存储回整数寄存器，以获得比较相等的元素的位图。 Popcount这些位图，并添加结果。

---

这让我想到了一个更好的主意：只需 3 次随机播放并按元素旋转即可完成所有比较：

{ 1d 1c 1b 1a }
  == == == ==   packed-compare with
{ 2d 2c 2b 2a }

{ 1a 1d 1c 1b }
  == == == ==   packed-compare with
{ 2d 2c 2b 2a }

{ 1b 1a 1d 1c }  # if dups didn't matter: do this shuffle on s2
  == == == ==   packed-compare with
{ 2d 2c 2b 2a }

{ 1c 1b 1a 1d } # if dups didn't matter: this result from { 1a ... }
  == == == ==   packed-compare with
{ 2d 2c 2b 2a }                                           { 2b ...

这只是 3 次随机播放，仍然进行了所有 16 次比较。诀窍是将它们与我们需要合并重复项的 OR 相结合，然后能够有效地计算它们。打包比较输出一个向量，其中每个元素 = 零或 -1（所有位设置），基于该位置的两个元素之间的比较。它旨在为 AND 或 XOR 创建一个有用的操作数来屏蔽一些向量元素，例如使 v1 += v2 & mask 以每个元素为条件。它也可以作为一个布尔真值。

通过将一个向量旋转 2 并将另一个向量旋转 1，然后在四个移位和未移位向量之间进行比较，所有 16 次比较只有 2 次随机播放是可能的。如果我们不需要消除重复，那就太好了，但既然我们这样做了，那么结果在哪里就很重要。我们不只是添加所有 16 个比较结果。

OR 将打包比较结果合并为一个向量。每个元素将根据 s2 的该元素是否在 s1 中有任何匹配来设置。 int _mm_movemask_ps (__m128 a) 将向量转换为位图，然后对位图进行计数。 （Nehalem or newer CPU required for popcnt，否则回退到具有 4 位查找表的版本。）

垂直 OR 处理 s1 中的重复项，但 s2 中的重复项是一个不太明显的扩展，需要更多的工作。我最终确实想到了一种速度不到两倍的方法（见下文）。

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

typedef struct four_points {
    int32_t a, b, c, d;
} four_points;
//typedef uint32_t four_points[4];

// small enough to inline, only 62B of x86 instructions (gcc 4.9.2)
static inline int match4_sse_noS2dup(const four_points *s1pointer, const four_points *s2pointer)
{
    __m128i s1 = _mm_loadu_si128((__m128i*)s1pointer);
    __m128i s2 = _mm_loadu_si128((__m128i*)s2pointer);
    __m128i s1b= _mm_shuffle_epi32(s1, _MM_SHUFFLE(0, 3, 2, 1));
    // no shuffle needed for first compare
    __m128i match = _mm_cmpeq_epi32(s1 , s2);  //{s1.d==s2.d?-1:0, 1c==2c, 1b==2b, 1a==2a }
    __m128i s1c= _mm_shuffle_epi32(s1, _MM_SHUFFLE(1, 0, 3, 2));
    s1b = _mm_cmpeq_epi32(s1b, s2);
    match = _mm_or_si128(match, s1b);  // merge dups by ORing instead of adding

    // note that we shuffle the original vector every time
    // multiple short dependency chains are better than one long one.
    __m128i s1d= _mm_shuffle_epi32(s1, _MM_SHUFFLE(2, 1, 0, 3));
    s1c = _mm_cmpeq_epi32(s1c, s2);
    match = _mm_or_si128(match, s1c);
    s1d = _mm_cmpeq_epi32(s1d, s2);

    match = _mm_or_si128(match, s1d);    // match = { s2.a in s1?,  s2.b in s1?, etc. }

    // turn the the high bit of each 32bit element into a bitmap of s2 elements that have matches anywhere in s1
    // use float movemask because integer movemask does 8bit elements.
    int matchmask = _mm_movemask_ps (_mm_castsi128_ps(match));

    return _mm_popcnt_u32(matchmask);  // or use a 4b lookup table for CPUs with SSE2 but not popcnt
}

查看在 s2 中消除重复的版本，以获取相同的代码，其中的行以更易读的顺序排列。我试图安排指令，以防 CPU 只是在执行指令之前几乎不解码指令，但是 gcc 将指令以相同的顺序排列，而不管你将内在函数放在什么顺序。

这非常快，如果在 128b 加载中没有存储转发停顿。如果您刚刚编写了具有四个 32 位存储的结构，则在接下来的几个时钟周期内运行此函数将在尝试以 128b 加载整个结构时产生停顿。见Agner Fog's site。如果调用代码已经在寄存器中包含 8 个值中的许多值，那么标量版本可能是一个胜利，即使对于只从内存中读取结构的微基准测试来说它会更慢。

我懒得为此进行循环计数，因为还没有完成重复处理。 IACA 表示，Haswell 可以以每 4.05 个时钟周期一次迭代的吞吐量和 17 个周期的延迟运行它（不确定这是否包括负载的内存延迟。有很多指令级并行性可用，并且所有指令都有单周期延迟，除了 movmsk(2) 和 popcnt(3))。如果没有 AVX，它会稍微慢一些，因为 gcc 选择了更差的指令顺序，并且仍然浪费了 movdqa 指令来复制向量寄存器。

使用 AVX2，这可以在 256b 个向量中并行执行两个 match4 操作。 AVX2 通常用作两个 128b 通道，而不是完整的 256b 向量。将您的代码设置为能够并行利用 2 个或 4 个 (AVX-512) match4 操作将为您在为这些 CPU 进行编译时带来收益。 s1s 或 s2s 不必连续存储，因此单个 32B 负载可以获得两个结构。 AVX2 具有相当快的加载 128b 到寄存器的上部通道。

---

处理s2 中的重复项

也许将 s2 与一个 shifted 而不是自身的旋转版本进行比较。

#### comparing S2 with itself to mask off duplicates
{  0 2d 2c 2b }
{ 2d 2c 2b 2a }     == == ==

{  0  0 2d 2c }
{ 2d 2c 2b 2a }        == ==

{  0  0  0 2d }
{ 2d 2c 2b 2a }           ==

嗯，如果零可以作为常规元素出现，我们可能还需要在比较之后进行字节移位，以将潜在的误报变成零。 如果在s1 中有一个无法出现的标记值，您可以移入其中的元素，而不是 0。（SSE 有 PALIGNR，它可以为您提供任何连续的 16B 窗口需要附加两个寄存器的内容。命名用于从两个对齐的负载模拟未对齐的负载的用例。因此您将拥有该元素的常量向量。）

---

更新：我想到了一个很好的技巧，可以避免对标识元素的需求。实际上，我们只需进行两次向量比较即可获得所有 6 次必要的 s2 与 s2 比较，然后将结果组合起来。

• 在两个向量的相同位置进行相同的比较可以让您将两个结果组合在一起，而无需在 OR 之前进行屏蔽。 （解决缺少标记值的问题）。

• 改组比较的输出，而不是 S2 的额外 shuffle&compare。这意味着我们可以在其他比较旁边完成d==a。

• 请注意，我们不仅限于将整个元素打乱。逐字节洗牌以从不同的比较结果中获取字节到单个向量元素中，并将 that 与零进行比较。 （这比我希望的要少，见下文）。

检查重复项会大大降低速度（尤其是在吞吐量方面，而不是在延迟方面）。所以你仍然最好在 s2 中安排一个永远不会匹配任何 s1 元素的哨兵值，你说这是可能的。我只介绍这个，因为我认为它很有趣。 （并为您提供一个选项，以防您有时需要不需要哨兵的版本。）

static inline
int match4_sse(const four_points *s1pointer, const four_points *s2pointer)
{
    // IACA_START
    __m128i s1 = _mm_loadu_si128((__m128i*)s1pointer);
    __m128i s2 = _mm_loadu_si128((__m128i*)s2pointer);
    // s1a = unshuffled = s1.a in the low element
    __m128i s1b= _mm_shuffle_epi32(s1, _MM_SHUFFLE(0, 3, 2, 1));
    __m128i s1c= _mm_shuffle_epi32(s1, _MM_SHUFFLE(1, 0, 3, 2));
    __m128i s1d= _mm_shuffle_epi32(s1, _MM_SHUFFLE(2, 1, 0, 3));

    __m128i match = _mm_cmpeq_epi32(s1 , s2);  //{s1.d==s2.d?-1:0, 1c==2c, 1b==2b, 1a==2a }
    s1b = _mm_cmpeq_epi32(s1b, s2);
    match = _mm_or_si128(match, s1b);  // merge dups by ORing instead of adding

    s1c = _mm_cmpeq_epi32(s1c, s2);
    match = _mm_or_si128(match, s1c);
    s1d = _mm_cmpeq_epi32(s1d, s2);
    match = _mm_or_si128(match, s1d);
    // match = { s2.a in s1?,  s2.b in s1?, etc. }

    // s1 vs s2 all done, now prepare a mask for it based on s2 dups

/*
 * d==b   c==a   b==a  d==a   #s2b
 * d==c   c==b   b==a  d==a   #s2c
 *    OR together -> s2bc
 *  d==abc     c==ba    b==a    0  pshufb(s2bc) (packed as zero or non-zero bytes within the each element)
 * !(d==abc) !(c==ba) !(b==a)  !0   pcmpeq setzero -> AND mask for s1_vs_s2 match
 */
    __m128i s2b = _mm_shuffle_epi32(s2, _MM_SHUFFLE(1, 0, 0, 3));
    __m128i s2c = _mm_shuffle_epi32(s2, _MM_SHUFFLE(2, 1, 0, 3));
    s2b = _mm_cmpeq_epi32(s2b, s2);
    s2c = _mm_cmpeq_epi32(s2c, s2);

    __m128i s2bc= _mm_or_si128(s2b, s2c);
    s2bc = _mm_shuffle_epi8(s2bc, _mm_set_epi8(-1,-1,0,12,  -1,-1,-1,8, -1,-1,-1,4,  -1,-1,-1,-1));
    __m128i dupmask = _mm_cmpeq_epi32(s2bc, _mm_setzero_si128());
    // see below for alternate insn sequences that can go here.

    match = _mm_and_si128(match, dupmask);
    // turn the the high bit of each 32bit element into a bitmap of s2 matches
    // use float movemask because integer movemask does 8bit elements.
    int matchmask = _mm_movemask_ps (_mm_castsi128_ps(match));

    int ret = _mm_popcnt_u32(matchmask);  // or use a 4b lookup table for CPUs with SSE2 but not popcnt
    // IACA_END
    return ret;
}

pshufb 需要 SSSE3。它和 pcmpeq（以及用于生成常量的 pxor）正在替换随机播放 (bslli(s2bc, 12))、OR 和 AND。

d==bc  c==ab b==a a==d = s2b|s2c
d==a   0     0    0    = byte-shift-left(s2b) = s2d0
d==abc c==ab b==a a==d = s2abc
d==abc c==ab b==a 0    = mask(s2abc).  Maybe use PBLENDW or MOVSS from s2d0 (which we know has zeros) to save loading a 16B mask.

__m128i s2abcd = _mm_or_si128(s2b, s2c);
//s2bc = _mm_shuffle_epi8(s2bc, _mm_set_epi8(-1,-1,0,12,  -1,-1,-1,8, -1,-1,-1,4,  -1,-1,-1,-1));
//__m128i dupmask = _mm_cmpeq_epi32(s2bc, _mm_setzero_si128());
__m128i s2d0 = _mm_bslli_si128(s2b, 12);  // d==a  0  0  0
s2abcd = _mm_or_si128(s2abcd, s2d0);
__m128i dupmask = _mm_blend_epi16(s2abcd, s2d0, 0 | (2 | 1));
//__m128i dupmask = _mm_and_si128(s2abcd, _mm_set_epi32(-1, -1, -1, 0));

match = _mm_andnot_si128(dupmask, match);  // ~dupmask & match;  first arg is the one that's inverted

我不能推荐MOVSS；它会在 AMD 上产生额外的延迟，因为它在 FP 域中运行。 PBLENDW 是 SSE4.1。 popcnt 在 AMD K10 上可用，但 PBLENDW 不可用（一些巴塞罗那核心的 PhenomII CPU 可能仍在使用中）。其实K10也没有PSHUFB，所以只需要SSE4.1和POPCNT，使用PBLENDW。 （或者使用 PSHUFB 版本，除非它会经常缓存丢失。）

避免从内存中加载向量常量的另一个选择是移动掩码 s2bc，并使用整数而不是向量操作。但是，它看起来会更慢，因为额外的移动掩码不是免费的，并且整数 ANDN 不可用。 BMI1 直到 Haswell 才出现，甚至 Skylake Celeron 和 Pentiums 也不会出现。 （Very annoying，IMO。这意味着 compilers can't start using BMI 更长。）

unsigned int dupmask = _mm_movemask_ps(cast(s2bc));
dupmask |= dupmask << 3;  // bit3 = d==abc.  garbage in bits 4-6, careful if using AVX2 to do two structs at once
        // only 2 instructions.  compiler can use lea r2, [r1*8] to copy and scale
dupmask &= ~1;  // clear the low bit

unsigned int matchmask = _mm_movemask_ps(cast(match));
matchmask &= ~dupmask;   // ANDN is in BMI1 (Haswell), so this will take 2 instructions
return _mm_popcnt_u32(matchmask);

AMD XOP 的 VPPERM（从两个源寄存器的任何元素中挑选字节）可以让字节洗牌取代合并 s2b 和 s2c 的 OR。

嗯，pshufb 并没有像我想象的那样节省我，因为它需要一个 pcmpeqd 和一个 pxor 来将寄存器归零。它还从内存中的常量加载其 shuffle 掩码，这可能会在 D-cache 中丢失。不过，这是我想出的最快的版本。

如果内联到循环中，可以使用相同的归零寄存器，从而节省一条指令。但是，OR 和 AND 可以在 port0（Intel CPU）上运行，它不能运行 shuffle 或 compare 指令。不过，PXOR 不使用任何执行端口（在英特尔 SnB 系列微架构上）。

我没有运行任何这些的真正基准，只有 IACA。

PBLENDW 和 PSHUFB 版本具有相同的延迟（22 个周期，为非 AVX 编译），但 PSHUFB 版本具有更好的吞吐量（每 7.1c 一个，而每 7.4c 一个，因为 PBLENDW 需要 shuffle 端口，并且已经有很多争论了。） IACA 表示，使用带有常量而不是 PBLENDW 的 PANDN 的版本也是每 7.4c 一个吞吐量，令人失望。 Port0 没有饱和，所以 IDK 为什么它和 PBLENDW 一样慢。

---

没有成功的旧想法。

为了让人们在将向量用于相关事物时寻找可以尝试的事物而受益。

使用向量对 s2 进行重复检查比检查 s2 与 s1 的工作量更大，因为如果使用向量进行比较，则与 4 次比较一样昂贵。 在比较之后需要改组或屏蔽，以在没有标记值的情况下消除误报，这很烦人。

目前的想法：

• 将s2 移动一个元素，并将其与自身进行比较。屏蔽误报以防止向 0 移动。将它们垂直或在一起，并使用它来 ANDN s1 与 s2 向量。

• 标量代码进行较少数量的 s2 与自身比较，构建一个位掩码以在 popcnt 之前使用。

• 广播s2.d 并对照s2（所有位置）进行检查。但这会将结果水平放置在一个向量中，而不是垂直放置在 3 个向量中。要使用它，也许PTEST / SETCC 为位图制作一个掩码（在popcount 之前应用）。 （PTEST 带有_mm_setr_epi32(0, -1, -1, -1) 的掩码，仅测试c,b,a，而不是d==d）。使用标量代码执行 (c==a | c==b) 和 b==a ，并将其组合成一个掩码。 Intel Haswell 及更高版本有 4 个 ALU 执行端口，但其中只有 3 个可以运行向量指令，因此混合中的一些标量代码可以填充端口 6。 AMD 在向量和整数执行资源之间有更多的分离。

• 随机播放s2 以某种方式完成所有必要的比较，然后随机播放输出。也许使用 movemask -> 4 位查找表？