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

Efficiently find least significant set bit in a large array?

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

有一条名为BSF 的汇编指令给出了最高设置位的位置(GCC 的__builtin_ctzll())。 所以在 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)

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

另请阅读:Find first set

您可以试试这个函数,您的编译器应该为您的 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;
}

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

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

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

相反,您想循环搜索非零向量,使用涉及 SSE4.1 ptest xmm0,xmm0 / [=17 的全向量检查=](3 微指令),或使用 SSE2 pcmpeqd v, zero / pmovmskb / cmp eax, 0xffff / je .loop(cmp/jcc 宏融合后 3 微指令)。 https://uops.info/

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

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

(如果预期的第一位位置不是非常高,那么是不值得的:如果第一个设置位在第一个向量中,在你加载的 2 个向量之间排序会更慢。5000 位只有 625 字节,或 19.5 个 AVX2 __m256i 向量。第一个设置位可能并不总是就在最后)

AVX2 版本:

这会检查成对的 32 字节向量(即整个缓存行)是否为非零值,如果找到则将其分类为一个 64 位位图以用于单个 CTZ 操作。额外的 shift/OR 会在关键路径中造成延迟,但希望我们能更快到达第一个 1 位。

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

(if 正文以 return 结尾,所以在 asm 中它实际上就像一个 if()break,或者实际上是一个 if()goto,因为它与未找到的 return -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) 操作的 CF 结果对全一向量进行检查,这样我们就可以检查所有元素是否匹配(即输入全为零)。不幸的是, - 不,我不认为我们可以在没有 vpor 的情况下有效地使用 vptest

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

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

没有 AVX

16 字节向量意味着我们不必处理 AVX2 洗牌的“通道内”行为。因此,我们可以结合 packssdwpacksswb 而不是 OR。包输入的高半部分中的任何设置位都会将结果符号饱和为 0x80 或 0x7f。 (所以带符号的饱和度是关键,而不是 unsigned packuswb 这会将带符号的负输入饱和到 0。)

但是,仅在 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 中 none 支持 AVX-512。即使是这样,您也可能希望避免使用 512 位向量,因为 SIMD instructions lowering CPU frequency,除非您的程序花费了 很多 的时间来执行此操作,并且您的数据在 L1d 缓存中很热,所以您可以真正受益,而不是仍然是 L2 缓存带宽的瓶颈。但即使使用 256 位向量,AVX-512 也有新的指令对此很有用:

  • 整数比较(vpcmpb/w/d/q) have a choice of predicate, so you can do not-equal instead of having to invert later with NOT. Or even test-into-register vptestmd 所以你不需要一个归零向量来比较。
  • compare-into-mask 有点像 pcmpeq + movmsk,除了结果在 k 寄存器中,还需要 kmovq rax, k0 才能 tzcnt.
  • kortest - 根据两个掩码寄存器的或非零设置FLAGS。所以搜索循环可以做 vpcmpd k0, ymm0, [rdi] / vpcmpd k1, ymm0, [rdi+32] / kortestw k0, k1

与多个输入数组

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

您可能希望在几个向量的块中执行此操作,乐观地希望早日在某处设置位。

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

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

L1d 缓存是 8 路关联的(在 Ice Lake 之前是 12 路),有限数量的 integer/pointer 寄存器可能使一次读取太多流成为一个坏主意。您可能不希望使编译器在指针内存中的实际数组上循环的间接级别。