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

Question

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

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

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

UPD：至于 CPU，我对此优化很感兴趣，以支持以下 CPU 阵容：
Intel Xeon E3-12XX、Intel Xeon E5-22XX/26XX/E56XX、Intel Core i3-5XX/4XXX/8XXX、Intel Core i5-7XX、Intel Celeron G18XX/G49XX（Intel Atom N2600、Intel Celeron N2807、Cortex-可选A53/72)

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

另请阅读：Find first set

Answer 1

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

相反，您希望遍历数组以搜索非零向量，使用涉及 SSE4.1 ptest xmm0,xmm0 / jz .loop (3 uops) 的全向量检查，或使用SSE2 pcmpeqd v, zero / pmovmskb / cmp eax, 0xffff / je .loop（cmp/jcc 宏融合后 3 微秒）。 https://uops.info/

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

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

（如果预期的第一位位置不是非常高，则不值得过于积极：如果第一个设置位在第一个向量，在你加载的 2 个向量之间排序会更慢。5000 位只有 625 个字节，或 19.5 个 AVX2 __m256i 向量。第一个设置位可能并不总是在最后）

AVX2 版本：

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

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

（if 主体以 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;
}

On Godbolt with 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，而不需要任何额外的前端微指令，只需要后端。只要编译器继续使用指针增量，而不是indexed addressing modes that would un-laminate。这将减少分支之后所需的工作量（这可能是错误的预测）。

要仍然使用vptest，您可能必须利用CF = (~dst & src == 0) 操作对全1 向量的CF 结果，因此我们可以检查所有元素是否匹配（即输入全为零） .不幸的是，Can PTEST be used to test if two registers are both zero or some other condition? - 不，我认为如果没有vpor，我们就无法有效地使用vptest。

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

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

没有 AVX

16 字节向量意味着我们不必处理 AVX2 shuffle 的“in-lane”行为。因此，我们可以使用packssdw 或packsswb 来代替OR。包输入的高半部分中的任何设置位将使结果符号饱和为 0x80 或 0x7f。 （所以有符号饱和度是关键，而不是 unsigned packuswb 它将使有符号负输入饱和为 0。）

但是，shuffle 仅在 Intel CPU 的端口 5 上运行，因此请注意吞吐量限制。例如，Skylake 上的 ptest 是 2 微指令，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));。即bsf eax, [rdi + rdx*2]。

使用 AVX-512：

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

• 整数比较 (vpcmpb/w/d/q) 可以选择谓词，因此您可以不等于，而不必稍后用 NOT 反转。甚至可以测试注册vptestmd，这样您就不需要一个归零向量来进行比较。
• compare-into-mask 有点像 pcmpeq + movmsk，除了结果在k 寄存器中，在tzcnt 之前仍然需要kmovq rax, k0。
• kortest - 根据两个非零掩码寄存器的 OR 设置 FLAGS。所以搜索循环可以做vpcmpd k0, ymm0, [rdi]/vpcmpd k1, ymm0, [rdi+32]/kortestw k0, k1

---

ANDing 多个输入数组

你提到你真正的问题是你有多达 20 个位数组，你想用 AND 与它们相交并在相交中找到第一个设置的位。

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

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

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

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

Answer 2

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

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

Answer 3

你可以试试这个功能，你的编译器应该为你的 CPU 优化这个代码。它不是超级完美，但它应该相对较快并且大部分是便携的。

PS length 应该能被 8 整除以达到最大速度

#include <stdio.h>
#include <stdint.h>

/* Returns the index position of the most significant bit; starting with index 0. */
/* Return value is between 0 and 64 times length. */
/* When return value is exact 64 times length, no significant bit was found, aka bf is 0. */
uint32_t offset_fsb(const uint64_t *bf, const register uint16_t length){
    register uint16_t i = 0;
    uint16_t remainder = length % 8;

    switch(remainder){
        case 0 : /* 512bit compare */
            while(i < length){
                if(bf[i] | bf[i+1] | bf[i+2] | bf[i+3] | bf[i+4] | bf[i+5] | bf[i+6] | bf[i+7]) break;
                i += 8;
            }
            /* fall through */

        case 4 : /* 256bit compare */
            while(i < length){
                if(bf[i] | bf[i+1] | bf[i+2] | bf[i+3]) break;
                i += 4;
            }
            /* fall through */

        case 6 : /* 128bit compare */    
            /* fall through */
        case 2 : /* 128bit compare */
            while(i < length){
                if(bf[i] | bf[i+1]) break;
                i += 2;
            }
            /* fall through */

        default : /* 64bit compare */
            while(i < length){
                if(bf[i]) break;
                i++;
            }
    }

    register uint32_t offset_fsb = i * 64;

    /* Check the last uint64_t if the last uint64_t is not 0. */
    if(bf[i]){
        register uint64_t s = bf[i];
        offset_fsb += 63;
        while(s >>= 1) offset_fsb--;
    }

    return offset_fsb;
}

int main(int argc, char *argv[]){
    uint64_t test[16];
    test[0] = 0;
    test[1] = 0;
    test[2] = 0;
    test[3] = 0;
    test[4] = 0;
    test[5] = 0;
    test[6] = 0;
    test[7] = 0;
    test[8] = 0;
    test[9] = 0;
    test[10] = 0;
    test[11] = 0;
    test[12] = 0;
    test[13] = 0;
    test[14] = 0;
    test[15] = 1;

    printf("offset_fsb = %d\n", offset_fsb(test, 16));

    return 0;
}