在大数组中有效地找到最低有效设置位？

Question

我有一个巨大的内存块（位向量），在一个内存页中大小为N位，考虑N平均为 5000，即 5k 位来存储一些标志信息。
在某个时间点（超级频繁 - 关键），我需要在整个大位向量中找到第一个位集。现在我每 64 个字都这样做，即在 ) 的帮助下__builtin_ctzll。但是当N增长并且搜索算法无法改进时，可以通过扩展内存访问宽度来扩展此搜索。这是几句话的主要问题

有一条被调用的汇编指令BSF 给出了最高设置位（GCC's __builtin_ctzll()）的位置。因此，在x86-64 arch 中，我可以在 64 位字中廉价地找到最高位。

但是通过内存宽度进行缩放呢？
例如，有没有办法用 128 / 256 / 512 位寄存器有效地做到这一点？
基本上我对一些 C API 函数来实现这个感兴趣，但也想知道这个方法是基于什么的。

UPD：至于 CPU，我对这种优化感兴趣，以支持以下 CPU 阵容：
英特尔至强 E3-12XX、英特尔至强 E5-22XX/26XX/E56XX、英特尔酷睿 i3-5XX/4XXX/8XXX、英特尔酷睿 i5- 7XX、英特尔赛扬 G18XX/G49XX（英特尔凌动 N2600、英特尔赛扬 N2807、Cortex-A53/72 可选）

PS在最终位扫描之前提到的算法中，我需要将k（平均 20-40）个N位向量与 CPU AND相加（AND 结果只是位扫描的准备阶段）。这也适用于内存宽度缩放（即比每 64 位字 AND 更有效）

另请阅读：查找第一组

Answer 1

这个答案是不同的，但如果你事先知道你将维护一个 B 位的集合，并且需要能够有效地设置和清除位，同时还要弄清楚哪个位是第一个设置的位，您可能想要使用像van Emde Boas 树或y-fast trie这样的数据结构。这些数据结构旨在存储小范围内的整数，因此您可以添加或删除要设置/清除的位的索引，而不是设置或清除单个位。它们非常快 - 您可以在 O(log log B) 时间内添加或删除项目，并且它们可以让您在 O(1) 时间内找到最小的项目。如图，如果 B ？50000，那么log log B大约是4。

我知道这并没有直接解决如何在巨大的位向量中找到最高位的问题。如果您的设置必须使用位向量，则其他答案可能会更有帮助。但是，如果您可以选择以不涉及位向量搜索的方式重新构建问题，那么这些其他数据结构可能更适合。

Answer 2

在整个向量中查找第一个设置位的最佳方法 (AFAIK) 包括查找第一个非零 SIMD 元素（例如字节或双字），然后对其使用位扫描。( __builtin_ctz/// bsf- tzcnt1 ffs) 。因此， ctz(vector) 本身并不是搜索数组的有用构建块，仅适用于循环之后。

相反，您希望使用涉及 SSE4.1 / (3 uops) 或 SSE2 // / （cmp/jcc 宏融合后的 3 uops）的 整体向量检查来循环搜索非零向量。https://uops.info/ptest xmm0,xmm0jz .looppcmpeqd v, zeropmovmskbcmp eax, 0xffffje .loop

一旦找到非零向量，pcmpeqb//movmskps就可以bsf 在其上找到双字索引，然后加载该双字及其bsf。将起始位位置 ( CHAR_BIT*4*dword_idx) 添加到bsf该元素内的位位置。这是一个相当长的延迟依赖链，包括整数 L1d 加载延迟。但由于您刚刚加载了向量，因此至少您可以相当有信心当您再次使用整数加载它时会在缓存中命中。（如果向量是动态生成的，那么最好还是存储/重新加载它并让存储转发工作，而不是尝试为vpermilps/movd或 SSSE3 pshufb//生成随机播放控制。movd）movzx ecx, al

strlen循环问题与or非常相似memchr，只不过我们拒绝单个值 (0) 并寻找其他值。尽管如此，我们仍然可以从手工优化的 asm strlen / memchr 实现（例如 glibc 的实现）中获得灵感，例如加载多个向量并进行一次检查以查看其中是否有任何向量具有所需的内容。（对于 strlen，pminub如果任何元素为 0，则与 组合得到 0。对于pcmpeqb比较结果，OR 对于 memchr）。出于我们的目的，我们想要的归约运算是 OR - 任何非零输入都会使输出非零，并且按位布尔运算可以在任何向量 ALU 端口上运行。

（如果预期的第一位位置不是很高，则不值得对此过于激进：如果第一个设置位位于第一个向量中，则在已加载的 2 个向量之间进行排序会更慢。5000位只有 625 个字节，或 19.5 个 AVX2__m256i向量。并且第一个设置位可能并不总是在末尾）

AVX2版本：

这会检查成对的 32 字节向量（即整个缓存行）是否非零，如果找到，则将其分类为单个 CTZ 操作的一个 64 位位图。额外的移位/或会导致关键路径中的延迟，但希望我们能早点到达第一个位。

使用 OR 将 2 个向量合并为 1 意味着了解 OR 结果的哪个元素非零并不是非常有用。我们基本上重做 if 内部的工作。这就是我们为保持实际搜索部分的微指令数量较低而付出的代价。

（if主体以 a 结尾return，因此在 asm 中，它实际上就像一个if()break，或者实际上if()goto是一个循环外，因为它到达的位置与未找到的返回 -1 因退出循环而到达不同的位置。）

// untested, especially the pointer end condition, but compiles to asm that looks good
// Assumes len is a multiple of 64 bytes

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

// aliasing-safe: p can point to any C data type
int bitscan_avx2(const char *p, size_t len /* in bytes */)
{
    //assert(len % 64 == 0);
    //optimal if p is 64-byte aligned, so we're checking single cache-lines
    const char *p_init = p;
    const char *endp = p + len - 64;
    do {
        __m256i v1 = _mm256_loadu_si256((const __m256i*)p);
        __m256i v2 = _mm256_loadu_si256((const __m256i*)(p+32));
        __m256i or = _mm256_or_si256(v1,v2);
        if (!_mm256_testz_si256(or, or)){        // find the first non-zero cache line
            __m256i v1z = _mm256_cmpeq_epi32(v1, _mm256_setzero_si256());
            __m256i v2z = _mm256_cmpeq_epi32(v2, _mm256_setzero_si256());
            uint32_t zero_map = _mm256_movemask_ps(_mm256_castsi256_ps(v1z));
            zero_map |= _mm256_movemask_ps(_mm256_castsi256_ps(v2z)) << 8;

            unsigned idx = __builtin_ctz(~zero_map);  // Use ctzll for GCC, because GCC is dumb and won't optimize away a movsx
            uint32_t nonzero_chunk;
            memcpy(&nonzero_chunk, p+4*idx, sizeof(nonzero_chunk));  // aliasing / alignment-safe load

            return (p-p_init + 4*idx)*8 + __builtin_ctz(nonzero_chunk);
        }
        p += 64;
    }while(p < endp);
    return -1;
}

在 Godbolt 上使用 clang 12 -O3 -march=haswell：

bitscan_avx2:
        lea     rax, [rdi + rsi]
        add     rax, -64                 # endp
        xor     ecx, ecx
.LBB0_1:                                # =>This Inner Loop Header: Depth=1
        vmovdqu ymm1, ymmword ptr [rdi]  # do {
        vmovdqu ymm0, ymmword ptr [rdi + 32]
        vpor    ymm2, ymm0, ymm1
        vptest  ymm2, ymm2
        jne     .LBB0_2                       # if() goto out of the inner loop
        add     ecx, 512                      # bit-counter incremented in the loop, for (p-p_init) * 8
        add     rdi, 64
        cmp     rdi, rax
        jb      .LBB0_1                  # }while(p<endp)

        mov     eax, -1               # not-found return path
        vzeroupper
        ret

.LBB0_2:
        vpxor   xmm2, xmm2, xmm2
        vpcmpeqd        ymm1, ymm1, ymm2
        vmovmskps       eax, ymm1
        vpcmpeqd        ymm0, ymm0, ymm2
        vmovmskps       edx, ymm0
        shl     edx, 8
        or      edx, eax             # mov ah,dl  would be interesting, but compilers won't do it.
        not     edx                  # one_positions = ~zero_positions
        xor     eax, eax                # break false dependency
        tzcnt   eax, edx             # dword_idx
        xor     edx, edx
        tzcnt   edx, dword ptr [rdi + 4*rax]   # p[dword_idx]
        shl     eax, 5               # dword_idx * 4 * CHAR_BIT
        add     eax, edx
        add     eax, ecx
        vzeroupper
        ret

这对于所有 CPU 来说可能不是最佳的，例如，也许我们可以为至少一个输入使用内存源vpcmpeqd，而不需要任何额外的前端微指令，只需要后端。只要编译器继续使用指针增量，而不是不会分层的索引寻址模式。这将减少分支之后所需的工作量（这可能会错误预测）。

要仍然使用vptest，您可能必须利用CF = (~dst & src == 0)针对全 1 向量进行操作的 CF 结果，因此我们可以检查所有元素是否匹配（即输入全为零）。不幸的是，PTEST 可以用来测试两个寄存器是否都为零或其他条件吗？- 不，我认为如果vptest没有vpor.

Clang 决定不在循环后实际减去指针，而是在搜索循环中做更多工作。:/ 循环为 9 uops（在cmp/宏融合之后jb），因此不幸的是，它每 2 个周期只能运行少于 1 次迭代。因此它只管理不到一半的 L1d 缓存带宽。

但显然单个数组不是你真正的问题。

不带AVX

16 字节向量意味着我们不必处理 AVX2 shuffle 的“车道内”行为。packssdw因此，我们可以用or来代替 OR 组合packsswb。包输入高半部分中的任何设置位都会将结果符号饱和为 0x80 或 0x7f。（因此有符号饱和是关键，而不是无符号packuswb，后者会使有符号负输入饱和为 0。）

但是，shuffle 仅在 Intel CPU 上的端口 5 上运行，因此请注意吞吐量限制。 ptest例如，在 Skylake 上是 2 uop、p5 和 p0，因此使用packsswb+ ptest+jz将限制每 2 个时钟进行一次迭代。但是pcmpeqd+pmovmskb不要。

不幸的是，在打包/组合之前pcmpeq单独使用每个输入会花费更多的微指令。但是会减少清理所需的工作量，并且如果循环退出通常涉及分支错误预测，则可能会减少总体延迟。

2x pcmpeqd=> packssdw=> pmovmskb=> not=>bsf会给你一个数字，你必须乘以 2 才能用作字节偏移量以获得非零双字。例如memcpy(&tmp_u32, p + (2*idx), sizeof(tmp_u32));。IE bsf eax, [rdi + rdx*2]。

使用 AVX-512：

您提到了 512 位向量，但您列出的 CPU 都不支持 AVX-512。即使是这样，您可能希望避免使用 512 位向量，因为SIMD 指令会降低 CPU 频率，除非您的程序花费大量时间来执行此操作，并且您的数据在 L1d 缓存中很热，因此您可以真正受益，而不是仍然在 L2 上成为瓶颈缓存带宽。但即使使用 256 位向量，AVX-512 也有对此有用的新指令：

• 整数比较 ( vpcmpb/w/d/q) 可以选择谓词，因此您可以不等于，而不必稍后使用 NOT 进行反转。或者甚至测试寄存器，vptestmd这样你就不需要一个归零向量来比较。

• Compare-into-mask 有点像 pcmpeq + movmsk，只不过结果在寄存器中k，仍然需要 akmovq rax, k0才能使用tzcnt。

• kortest- 根据两个掩码寄存器的或非零设置FLAGS。所以搜索循环可以做vpcmpd k0, ymm0, [rdi]//vpcmpd k1, ymm0, [rdi+32]kortestw k0, k1

• vplzcntd（或q） - 与 SIMD 结合isolate_lowest = v &= -v，可以找到最低设置位的位置（在每个 SIMD 向量中）。对于非零输入，bit_index = 31-lzcnt = 31^lzcnt。

• vpcompressq/ d- reg-reg 版本 ( https://uops.info ) 在 Intel 和 Zen 4 上为 2 uops。接下来vmovq eax, ymm0，这可以提取最低的非零元素（给定比较掩码），其延迟可能比tzcnt掩码上的标量要低，以索引另一个负载。

  但您仍然需要该标量tzcnt来找出要添加到双字内位索引中的内容，因此这会花费额外的微指令来缩短关键路径延迟。例如

// untested and worse for throughput, probably better for latency.
// Just writing it out to see what it looks like

// after finding a v  with a a non-zero bit somewhere:
  __mmask8 nzmask = _mm256_test_epi32_mask(v,v);  // true for non-zero elements
  __m256i bit_in_dword_lzcnt = _mm256_lzcnt_epi32(v & -v);  // lzcnt of the lowest set bit
  __m256i tmp = _mm256_maskz_compress_epi32(nzmask, bit_in_dword_lzcnt);  // low element has the lzcnt we want

  unsigned bit_idx = _tzcnt_u32(nzmask)*32;
  bit_idx += 31^_mm_cvtsi128_si32(_mm256_castsi256_si128(tmp)); // vmovd + xor to do 31-lzcnt more cheaply.

根据uops.info，vpcompressd英特尔的延迟从掩码到输出是6个周期，但从向量输入到向量输出只有3个周期。所以我猜第一个 uop 只是将掩码预处理为vpermd随机播放控件。

在 Zen 4 上，对于 256 位向量宽度，从向量输入到输出需要 4 个周期，从掩码到输出需要 8 个周期。对于 512 位，8:9。

来自矢量输入的时间vplzcntd(v & -v)比vptestmd(v)获取掩码花费的时间更长，因此效果很好。

对多个输入数组进行 AND 运算

您提到您真正的问题是您有最多 20 个位数组，并且您想将它们与 AND 相交并找到交集中的第一个设置位。

您可能希望在几个向量的块中执行此操作，乐观地希望尽早在某个地方有一个集合位。

AND 组由 4 或 8 个输入组成，通过 OR 累加结果，这样您就可以判断每个输入可能有 4 个向量的块中是否有 1。（如果没有任何 1 位，则在仍然加载指针的同时再执行 4 个向量、64 或 128 字节的块，因为如果您现在转到其他输入，交集肯定是空的）。调整这些块大小取决于 1 的稀疏程度，例如，可能始终以 6 或 8 个向量的块工作。不过，2 的幂数字很好，因为您可以将分配填充到 64 或 128 字节的倍数，这样您就不必担心提前停止。）

（对于奇数个输入，可以将同一个指针两次传递给需要 4 个输入的函数，而不是为每个可能的数字分派到循环的特殊版本。）

L1d 缓存是 8 路关联的（在 Ice Lake 之前是 12 路），并且有限数量的整数/指针寄存器可能会使尝试一次读取太多流成为一个坏主意。您可能也不希望间接级别使编译器循环遍历指针内存中的实际数组。