在比特流中查找比特模式的最快方法

39

我需要在位流中扫描一个16位的字。不能保证对齐在字节或字边界上。

实现这个最快的方式是什么?有各种各样的暴力方法;使用表格和/或移位,但是否有任何"位操作快捷方式"可以通过为每个字节或字到达时给出包含标志的yes/no/maybe结果来减少计算次数?

C代码、指令、x86机器码均有意义。


有多少搜索模式?这些模式在编译时已知吗? - LukeH
这个词是固定的吗?也就是说,算法能够针对这个词进行定制吗?还是你需要一个通用算法,可以在所有65536个可能的单词中都表现得相当不错? - MSalters
它可以在编译时固定。 - user172783
并且在比特流中是唯一的。 - user172783
哇,我真的很喜欢这个问题的回答,继续来吧!!! - AntonioMO
13个回答

27

有时使用简单的暴力方法是好的。

我认为预先计算单词的所有移位值,并将它们放入16个整数中,这样你就得到了一个类似于以下数组(假设intshort宽两倍):

 unsigned short pattern = 1234;
 unsigned int preShifts[16];
 unsigned int masks[16];
 int i;
 for(i=0; i<16; i++)
 {
      preShifts[i] = (unsigned int)(pattern<<i);  //gets promoted to int
      masks[i] = (unsigned int) (0xffff<<i);
 }

然后,对于从流中获取的每个无符号短整型,将该短整型和前一个短整型组成一个整数,并将该无符号整数与16个无符号整数进行比较。如果其中任何一个匹配,则表示找到了一个。

基本上就是这样:

  int numMatch(unsigned short curWord, unsigned short prevWord)
  {
       int numHits = 0;
       int combinedWords = (prevWord<<16) + curWord;

       int i=0;
       for(i=0; i<16; i++)
       {
             if((combinedWords & masks[i]) == preShifsts[i]) numHits++;
       }
       return numHits;
  }

请注意,这可能意味着在相同的位上检测到多个模式时会产生多次命中:
例如,在32位0和您想要检测的16个0的模式下,这将意味着该模式被检测到16次!
如果这段代码能够按照既定方式编译,那么它每输入一个单词就需要进行16次检查。对于每个输入比特位,它会执行一次 "&" 操作和 "==" 操作,以及分支或其他条件递增操作。此外,还会为每个比特位查找掩码。
查找表是不必要的;通过对 "combined" 进行右移操作,我们可以获得更高效的汇编代码,如 在另一个答案 中所示,该答案还展示了如何在 x86 上使用 SIMD 进行向量化。

1
这更符合我的思路,并且很可能是最好的选择,特别是如果SIMD内在函数允许跨多个字进行多个比较。 "标志"在位流中是非零且唯一的。 - user172783
2
很棒。如果你使用SIMD,通过将模式在巨大的寄存器中“并排”几次,你可以潜在地更快地进行搜索。这样,您可以同时比较4个无符号整数之类的内容。 - Toad
2
迈克尔:我认为你没有完全理解逻辑。16位模式被移位16次进入16个无符号整数,从而创建了16位字的所有可能变化。通过每次获取16位并进行掩码处理,然后将其与16个可能的模式进行比较,您已经完成了所有需要的检查。 - Toad
我看到上述算法存在问题,如果在短字节边界(对齐)处出现模式匹配,则可能会导致(错误的)重复计数。例如,如果模式是xaabb,流是x1122aabb3344。第一次调用numMatch(xaabb,x1122)将返回1。第二次调用numMatch(x3344,xaabb)将返回1。总匹配次数为2。 - Ahmed A
由于我们通过将0位移动到包括15位来预先计算模式,我不明白您建议如何匹配16位进入流中的模式。 - Toad
显示剩余9条评论

17
这里有一个技巧,可以将搜索速度提高32倍,如果使用{0,1}字母表的Knuth-Morris-Pratt算法或reinier的方法都不够快。你可以先使用一个具有256个条目的表来检查你的比特流中的每个字节是否包含在你要查找的16位单词中。使用以下代码可以获得该表:
unsigned char table[256];
for (int i=0; i<256; i++)
  table[i] = 0; // initialize with false
for (i=0; i<8; i++)
  table[(word >> i) & 0xff] = 1; // mark contained bytes with true

使用该方法,您可以找到位流中匹配的可能位置:

for (i=0; i<length; i++) {
  if (table[bitstream[i]]) {
    // here comes the code which checks if there is really a match
  }
}

由于256个表项中最多只有8个不为零,因此平均而言,您只需要查看每32个位置的内容。只需对这个字节(与前面和后面的字节组合)使用reinier建议的位操作或掩码技术来检查是否匹配。

该代码假定您使用小端字节顺序。字节中的位顺序也可能是一个问题(对于已经实现CRC32校验和的人来说也是众所周知的问题)。


1
聪明的加速方式! - Toad
如果空间有限,该表可以缩小为32字节的位数组。 - Jeremy

10

我建议使用大小为256的3个查找表的解决方案。这对于大型位流非常有效。此解决方案需要取样3个字节进行比较。下图显示了16位数据在3个字节中的所有可能排列方式,每个字节区域以不同的颜色显示。

alt text http://img70.imageshack.us/img70/8711/80541519.jpg

现在,在第一个样本中将处理1到8的检查,在下一个样本中处理9到16,依此类推。现在,当我们搜索Pattern时,我们将找到此Pattern的所有8种可能的排列方式(如下),并将其存储在3个查找表(左、中和右)中。

初始化查找表:

让我们以0111011101110111作为要查找的Pattern的示例。现在考虑第4个排列。左侧部分将是XXX01110。用00010000填充由左侧部分(XXX01110)指向的左侧查找表的所有行。 1表示输入Pattern的排列的起始位置。因此,左侧查找表的以下8行将填充为16(00010000)。

00001110
00101110
01001110
01101110
10001110
10101110
11001110
11101110

安排的中间部分将是11101110。在中间查找表中通过索引(238)指向的原始指针将被填充为16(00010000)。
现在,安排的右侧部分将是111XXXXX。所有索引为111XXXXX的原始数据(32个原始数据)将被填充为16(00010000)。
我们在填充时不应覆盖查找表中的元素。相反,执行按位OR操作以更新已经填充的原始数据。在上面的例子中,由第3个安排编写的所有原始数据将由第7个安排如下更新。

alt text

因此,左侧查找表中索引为XX011101,中间查找表中为11101110,右侧查找表中为111XXXXX的原始数据将通过第7次排列更新为00100010

搜索模式:

取三个字节的样本。按照以下方式找到计数,其中Left是左侧查找表,Middle是中间查找表,Right是右侧查找表。

Count = Left[Byte0] & Middle[Byte1] & Right[Byte2];

{{Count}}中的1的数量给出了在采样中匹配{{Pattern}}的数量。

我可以提供一些已测试的样本代码。

初始化查找表:

    for( RightShift = 0; RightShift < 8; RightShift++ )
    {
        LeftShift = 8 - RightShift;

        Starting = 128 >> RightShift;

        Byte = MSB >> RightShift;

        Count = 0xFF >> LeftShift;

        for( i = 0; i <= Count; i++ )
        {
            Index = ( i << LeftShift ) | Byte;

            Left[Index] |= Starting;
        }

        Byte = LSB << LeftShift;

        Count = 0xFF >> RightShift;

        for( i = 0; i <= Count; i++ )
        {
            Index = i | Byte;

            Right[Index] |= Starting;
        }

        Index = ( unsigned char )(( Pattern >> RightShift ) & 0xFF );

        Middle[Index] |= Starting;
    }

搜索模式:

Data 是流缓冲区,Left 是左侧查找表,Middle 是中间查找表,Right 是右侧查找表。

for( int Index = 1; Index < ( StreamLength - 1); Index++ )
{
    Count = Left[Data[Index - 1]] & Middle[Data[Index]] & Right[Data[Index + 1]];

    if( Count )
    {
        TotalCount += GetNumberOfOnes( Count );
    }
}

限制:

以上循环无法检测到流缓冲区末尾放置的模式。需要在循环后添加以下代码以克服此限制。

Count = Left[Data[StreamLength - 2]] & Middle[Data[StreamLength - 1]] & 128;

if( Count )
{
    TotalCount += GetNumberOfOnes( Count );
}

优点:

该算法只需要N-1个逻辑步骤就能在包含N个字节的数组中找到一个模式。唯一的开销是最初填充查找表,这在所有情况下都是恒定的。因此,对于搜索大量字节流非常有效。


2
太好了,解决方案仍在不断涌现。顺便说一下,你的名字听起来很棒。可能是《星球大战》电影中一个角色的名字;^) - Toad
2
我能看到两个优化点。第一个是在查找左右之前先检查中间,如果中间不好看,则跳过左右检查。第二个优化是使用一个256x32位表而不是三个256x8表。这将用一个长字延迟替换三个一字节的延迟(匹配公式为(byte2_lookup >> 16)&(byte1_lookup >> 8)&(byte0_lookup))。请注意,如果查找值为零,则可以跳过一个字节,前提是在找到非零查找值时准备回溯。 - supercat
“raw” 应该是 “row” 吗?术语看起来很奇怪,但查找表的基本思想是检查字节是否在某个偏移量上匹配。我们可以尝试将其适应于 4 位块的 SIMD 查找,例如使用 4 位索引从字节向量中并行执行 16 次查找的 SSSE3 pshufb - Peter Cordes

9

4
“两个字符的字母表”意味着将字节流转换为比特流。 - MSalters
2
Knuth自动将分数加到答案中。;^) 我认为他们的算法与位匹配无关。 - Toad
6
KMP算法的性能表现不佳,因为它依赖于跳过字符。然而,在搜索10……时,你无法跳过任何字符。当遇到11……时,你将遇到不匹配的情况,但是字符串可能是110……,所以必须将搜索串中的第二个1与模式的第一个1进行比较。如果你的模式是11111……,KMP算法可以很好地工作;每当遇到0,你可以跳过5个字符。KMP算法在字符集较大的情况下运行最佳,例如Unicode。对于只有两个字符的字符集来说,KMP算法是最差的情况(但其他一些搜索算法也会因字符集较小而表现不佳)。 - MSalters
5
对于他遇到的问题而言,使用KMP算法有些过度了,因为KMP算法适用于任意长度的单词,而他的单词长度固定为16位。针对这个问题来说,KMP算法过于灵活,不如采用更简单的方法更高效。 - A. Levy
一次只提取一个位的运行时开销将与像Reinier这样的解决方案中的整个工作量相当,而此时您甚至还没有开始KMP工作。除此之外,仅有两个字符和短的固定长度搜索模式的字母表是KMP的最坏情况。请记住,生活不仅仅是大O分析。 - Alan
显示剩余2条评论

7
我会实现一个拥有16个状态的状态机。
每个状态代表符合模式的接收位数。如果下一个接收的比特与模式的下一个比特相符,则机器步进到下一个状态。如果不是这种情况,则机器回到第一个状态(或者如果模式的开头可以匹配较少的接收位数,则回到另一个状态)。
当机器达到最后一个状态时,这表示已在比特流中识别出该模式。

这将与输入比特流一样快。 - mouviciel
对于任何实际应用来说,这种方法的速度虽然可能不如那些尝试同时匹配多个版本模式的版本快,但已经足够快了。 - Mark Bessey
1
这可能会非常快,而且还有另一个优点,就是非常简单,因此边界条件错误的可能性更小。唯一棘手的部分是处理部分匹配-您必须将状态“重新设置回流的某个位置”。简单的方法只是备份数据流,可能会导致在病态流(或模式)中变慢,其中您会获得许多15位匹配项。可以在状态机中构建正确的状态以重新启动匹配,而无需重新检查位,但这将使状态机变得更加棘手(我认为)。 - Michael Burr
在一般情况下,这很棘手。但如果模式是问题的常量,状态重置可以硬编码。或者您可以将棘手的部分放在初始化步骤中,并使用以模式为参数的状态机生成函数。 - mouviciel
@ziggystar:你能详细说明一下吗? - mouviciel
@mouviciel 对不起,我没有仔细阅读你的答案。我已经删除了我的错误评论。 - ziggystar

4
一种更简单的实现@Toad's simple brute-force algorithm that checks every bit-position的方法是将数据移位到位,而不是移动掩码。在循环内部只需右移combined >>= 1并比较低16位即可,无需任何数组。 (可以使用固定掩码或转换为uint16_t)。
(在多个问题中,我注意到创建掩码往往不如仅移出不需要的位效率高。)
(正确处理uint16_t数组的最后一个16位块,特别是奇数大小的字节数组的最后一个字节,留给读者作为练习。)
// simple brute-force scalar version, checks every bit position 1 at a time.
long bitstream_search_rshift(uint8_t *buf, size_t len, unsigned short pattern)
{
        uint16_t *bufshort = (uint16_t*)buf;  // maybe unsafe type punning
        len /= 2;
        for (size_t i = 0 ; i<len-1 ; i++) {
                //unsigned short curWord = bufshort[i];
                //unsigned short prevWord = bufshort[i+1];
                //int combinedWords = (prevWord<<16) + curWord;

                uint32_t combined;                                // assumes little-endian
                memcpy(&combined, bufshort+i, sizeof(combined));  // safe unaligned load

                for(int bitpos=0; bitpos<16; bitpos++) {
                        if( (combined&0xFFFF) == pattern)     // compiles more efficiently on e.g. old ARM32 without UBFX than (uint16_t)combined
                                return i*16 + bitpos;
                        combined >>= 1;
                }
        }
        return -1;
}

这比在大多数ISA(如x86、AArch64和ARM)上从数组加载掩码在最近的gcc和clang中编译效率显著提高。

编译器将循环完全展开16次,以便使用带有立即操作数的位字段提取指令(例如ARM的ubfx无符号位字段提取或PowerPC rwlinm左旋+立即掩码一个位范围),将16位提取到32或64位寄存器的底部,在那里可以进行常规的比较和分支。实际上没有1个右移的依赖链。

在x86上,CPU可以执行忽略高位的16位比较,例如在右移edx中的combined之后执行cmp cx,dx

一些ISA的编译器能够像@Toad版本那样很好地处理,例如PowerPC的clang能够优化掉掩码数组,使用rlwinm使用立即数来屏蔽combined的16位范围,并将所有16个预移位模式值保留在16个寄存器中,因此无论rlwinm是否具有非零旋转计数,都只是rlwinm / compare / branch。 但右移版本不需要设置16个tmp寄存器。 https://godbolt.org/z/8mUaDI


AVX2暴力破解

有(至少)两种方法可以实现:

  • 广播一个单个dword并使用可变移位来检查其所有位位置,然后再继续。可能很容易确定您找到匹配的位置。(如果您想要计算所有匹配项,则可能不太好。)
  • 向量加载,并同时迭代多个窗口数据的位位置。 可能使用不对齐的负载开始于相邻的字(16位),以获取dword(32位)窗口。否则,您将不得不跨越128位车道进行洗牌,最好具有16位粒度,并且这将需要2个指令而没有AVX512。
使用 64 位元素移位而不是 32 位,我们可以检查多个相邻的 16 位窗口,而不总是忽略上部的 16 位(其中会移入零)。但我们仍然在 SIMD 元素边界处中断,在那里会移入零,而不是来自更高地址的实际数据。(未来解决方案:AVX512VBMI2 双重移位像 VPSHRDW,这是 SHRD 的 SIMD 版本。)也许无论如何都值得这样做,然后回来处理每个 64 位元素顶部遗漏的 4 个 16 位元素,在一个 __m256i 中合并多个向量的剩余部分。
// simple brute force, broadcast 32 bits and then search for a 16-bit match at bit offset 0..15

#ifdef __AVX2__
#include <immintrin.h>
long bitstream_search_avx2(uint8_t *buf, size_t len, unsigned short pattern)
{
    __m256i vpat = _mm256_set1_epi32(pattern);

    len /= 2;
    uint16_t *bufshort = (uint16_t*)buf;
    for (size_t i = 0 ; i<len-1 ; i++) {
        uint32_t combined; // assumes little-endian
        memcpy(&combined, bufshort+i, sizeof(combined));  // safe unaligned load
        __m256i v = _mm256_set1_epi32(combined);
//      __m256i vlo = _mm256_srlv_epi32(v, _mm256_set_epi32(7,6,5,4,3,2,1,0));
//      __m256i vhi = _mm256_srli_epi32(vlo, 8);

        // shift counts set up to match lane ordering for vpacksswb

        // SRLVD cost: Skylake: as fast as other shifts: 1 uop, 2-per-clock
        // * Haswell: 3 uops
        // * Ryzen: 1 uop, but 3c latency and 2c throughput.  Or 4c / 4c for ymm 2 uop version
        // * Excavator: latency worse than PSRLD xmm, imm8 by 1c, same throughput. XMM: 3c latency / 1c tput.  YMM: 3c latency / 2c tput.  (http://users.atw.hu/instlatx64/AuthenticAMD0660F51_K15_BristolRidge_InstLatX64.txt)  Agner's numbers are different.
        __m256i vlo = _mm256_srlv_epi32(v, _mm256_set_epi32(11,10,9,8,    3,2,1,0));
        __m256i vhi = _mm256_srlv_epi32(v, _mm256_set_epi32(15,14,13,12,  7,6,5,4));

        __m256i cmplo = _mm256_cmpeq_epi16(vlo, vpat);  // low 16 of every 32-bit element = useful
        __m256i cmphi = _mm256_cmpeq_epi16(vhi, vpat);

        __m256i cmp_packed = _mm256_packs_epi16(cmplo, cmphi); // 8-bit elements, preserves sign bit
        unsigned cmpmask = _mm256_movemask_epi8(cmp_packed);
        cmpmask &= 0x55555555;  // discard odd bits
        if (cmpmask) {
            return i*16 + __builtin_ctz(cmpmask)/2;
        }
    }
    return -1;
}
#endif

这对于通常很快找到命中的搜索非常有用,特别是在数据的前32个字节内。对于大型搜索来说也不错(但仍然是纯暴力,每次只检查一个单词),在Skylake上,同时检查多个窗口的16个偏移量可能不比它更差。
这是针对Skylake进行调整的,在其他CPU上,其中可变移位效率较低,您可以考虑仅对0..7偏移量进行1次可变移位,然后通过移动来创建8..15偏移量。或者完全使用其他方法。

This compiles surprisingly well with gcc/clang (on Godbolt), with an inner loop that broadcasts straight from memory. (Optimizing the memcpy unaligned load and the set1() into a single vpbroadcastd)

此外,Godbolt链接中还包括一个测试main,它在一个小数组上运行。 (自上次微调以来我可能没有测试过,但我之前测试过,打包+位扫描功能确实有效。)
## clang8.0  -O3 -march=skylake  inner loop
.LBB0_2:                                # =>This Inner Loop Header: Depth=1
        vpbroadcastd    ymm3, dword ptr [rdi + 2*rdx]   # broadcast load
        vpsrlvd ymm4, ymm3, ymm1
        vpsrlvd ymm3, ymm3, ymm2             # shift 2 ways
        vpcmpeqw        ymm4, ymm4, ymm0
        vpcmpeqw        ymm3, ymm3, ymm0     # compare those results
        vpacksswb       ymm3, ymm4, ymm3     # pack to 8-bit elements
        vpmovmskb       ecx, ymm3            # scalar bitmask
        and     ecx, 1431655765              # see if any even elements matched
        jne     .LBB0_4            # break out of the loop on found, going to a tzcnt / ... epilogue
        add     rdx, 1
        add     r8, 16           # stupid compiler, calculate this with a multiply on a hit.
        cmp     rdx, rsi
        jb      .LBB0_2                    # } while(i<len-1);
        # fall through to not-found.

这是8个工作uops加上3个循环开销uops(假设在Haswell/Skylake上使用and/jne和cmp/jb的宏融合)。在AMD上,256位指令需要多个uops,因此会更多。


或者当然可以使用纯右移立即数将所有元素向右移动1位,并且同时检查多个窗口而不是在同一个窗口中使用多个偏移量。

如果没有高效的变量移位(特别是根本没有AVX2),对于大型搜索来说,即使需要更多的工作来确定第一个命中的位置,这也比使用多个偏移量更好。 (在找到除最低元素以外的其他地方的命中后,您需要检查所有较早窗口的所有剩余偏移量。)


看起来加速效果很不错。点赞实际检查汇编代码。我总是听到声称“编译器很聪明,会做一些聪明的事情”,但过去我已经看到这并非总是如此。很高兴看到当前的编译器实际上相当聪明。 - Toad
@Toad:是的,编译器通常会生成不那么糟糕的代码,但有些人对编译器的信任程度太过离谱了。它们经常需要一些手把手的指导!是的,clang优化掉了一些数组,这是比较聪明的做法。但并非所有架构都能实现这一点。对于一些架构,如x86,它们会自动向量化数组初始化,然后在标量循环中使用它。即使ICC也无法自动向量化搜索本身。(而且ICC的标量循环远非聪明;除了自动向量化之外,并不出色)。因此,AVX2版本必须手动进行向量化处理。 - Peter Cordes

4

atomice的

看起来很好,直到我考虑到Luke和MSalter关于特定信息的更多要求。

事实证明,特定情况可能比KMP更快。 KMP文章链接到

对于搜索模式为“ AAAAAA”的特定情况,可以使用。 对于多模式搜索,

可能最合适。

您可以在此处找到更多入门讨论


2
仍然...这些都涉及文本搜索,而没有涉及位级别的内容。将比特流先转换为字节流可能非常昂贵(内存/处理器时间)和不切实际。 - Toad
Reinier,我在这里专注于字符串搜索算法。你指出按位掩码操作并不是免费的。目前,我假设它们对于每个算法来说是相对昂贵的。然而,我的主要观点是应用程序的具体情况可能使我们击败KMP算法。 - Ewan Todd
KMP 只能一次测试一个字母。使用位运算,您可以一次测试 16 个字母(或者使用 SIMD 64)。这将使 KMP 或任何其他基于字母的算法变得无用。 - Toad
1
我不明白你坚持认为KMP仅限于字节大小操作的理由。我想到了一种按位运算的KMP算法,它在不匹配时向右移动一位。 - Ewan Todd
2
请搜索“将Knuth-Morris-Pratt算法应用于Huffman编码文本的模式匹配适应性”。 - Ewan Todd

3

对于通用的、非SIMD算法,你很难比这个更好:

unsigned int const pattern = pattern to search for
unsigned int accumulator = first three input bytes

do
{
  bool const found = ( ((accumulator   ) & ((1<<16)-1)) == pattern )
                   | ( ((accumulator>>1) & ((1<<16)-1)) == pattern );
                   | ( ((accumulator>>2) & ((1<<16)-1)) == pattern );
                   | ( ((accumulator>>3) & ((1<<16)-1)) == pattern );
                   | ( ((accumulator>>4) & ((1<<16)-1)) == pattern );
                   | ( ((accumulator>>5) & ((1<<16)-1)) == pattern );
                   | ( ((accumulator>>6) & ((1<<16)-1)) == pattern );
                   | ( ((accumulator>>7) & ((1<<16)-1)) == pattern );
  if( found ) { /* pattern found */ }
  accumulator >>= 8;

  unsigned int const data = next input byte
  accumulator |= (data<<8);
} while( there is input data left );

你的累加器在初始化时需要3个输入字节 ;^) - Toad

3
您可以使用快速傅里叶变换来处理极大的输入(n的值),以在O(n log n)的时间内查找任何位模式。计算位掩码与输入的交叉相关性。序列x和掩码y的交叉相关性,分别具有大小为n和n',定义如下:
R(m) = sum  _ k = 0 ^ n' x_{k+m} y_k

然后,如果您的位模式与掩码完全匹配并且 R(m) = Y,其中 Y 是您的位掩码中 1 的总数,则会发生匹配。

因此,如果您正在尝试匹配位模式

[0 0 1 0 1 0]

[ 1 1 0 0 1 0 1 0 0 0 1 0 1 0 1]

如果您想使用口罩,则必须戴上它。

[-1 -1  1 -1  1 -1]

掩码中的-1保证这些位置必须为0。

您可以使用FFT在O(n log n)时间内实现交叉相关。

我认为KMP具有O(n + k)的运行时间,因此它比这个更好。


2
如果这个问题可以在O(n)时间内解决,你会使用傅里叶吗? - Luka Rahne

3
似乎SIMD指令是一个很好的运用。SSE2为同时计算多个整数添加了一些整数指令,但我无法想象有多少解决方案不涉及大量位移,因为您的数据不会被对齐。这实际上听起来像是FPGA应该做的事情。

1
在软件中模拟FPGA的工作可能是一个有趣的起点。 - user172783
1
所有FPGA所做的就是移位和比较,电路并不关心字节对齐,它可以在比较当前位时同时移入下一位。 - ajs410

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接