[翻译] 优化x86上小规模负载的CRC32

这位作者是游戏机硬件爱好者，业余时间从事游戏机模拟器的开发，原文是针对小规模的CRC32的优化内容

---

最近在一组夜班中，我疲惫的头脑想了解如何通过无进位乘法运算来实现CRC32。这似乎是出于以下几种想法:

1. x86支持pclmulqdq指令，实现64bit * 64bit -> 127bit的无进位乘法
2. 尽管x86有针对Castagnoli多项式0x1EDC6F41的SSE4.2的crc32c指令，但是它好像没有针对其他多项式的类似的指令
3. AArch64系统提供了一个使用ISO/ANSI/gzip/PNG多项式0x04C11DB7的CRC32指令
4. 现在的AArch64模拟器回退到缓慢的软件实现当遇到了上述提到的指令

我对CRC32有一个模糊的理解，所以我知道这是这是可能的，但是当时的理解不足以使实现变得显而易见。回想起来，该解决方案应该是显而易见的，但是我还是在为了后人记录一下。

需求实施

主要动机是重新实现AArch64架构下的crc32指令集，因此让我们简单看下它们。它们是一系列针对小负载的指令:

CRC32B <Wd>, <Wn>, <Wm>
CRC32H <Wd>, <Wn>, <Wm>
CRC32W <Wd>, <Wn>, <Wm>
CRC32X <Wd>, <Wn>, <Xm>

针对所有的这些指令，n寄存器是累计的crc32值，m寄存器是当前的输入值。可以在AArch manual查看这些指令的ASL实现:

bits(32)    acc  = X[n]; // accumulator
bits(size)  val  = X[m]; // input value
bits(32)    poly = 0x04C11DB7;

bits(32+size) tempacc = BitReverse(acc):Zeros(size);
bits(size+32) tempval = BitReverse(val):Zeros(32);

// Poly32Mod2 on a bitstring does a polynomial Modulus over {0, 1} operation
X[d] = BitReverse(Poly32Mod2(tempacc EOR tempval, poly));

从上面我们能看出我们需要实现8,16,32和64位负载的crc32。

循环冗余校验(CRC)运算

\(crc(m(x)) = x ^ {degree(p(x))} * m(x) \pmod {p(x)}\)

如上所示，一个消息的CRC是附带了0的消息，表示为生成多项式的模的多项式。附加到消息的0的数量是多项式的度。

多项式的系数是Galois field GF(2)的元素。 所以，系数的加法等同于xor，并且两个多项式的乘法等同于无进位的乘法。两个多项式的除法等于无进位的除法(模2除法)。

对于ISO CRC32而言，\(degree(p(x))=32\)，并且生成多项式为:

\[p(x) = x ^ {32} + x ^ {26} + x ^ {23} + x ^ {22} + x ^ {16} + x ^ {12} + x ^ {11} + x ^ {10} + x ^ 8 + x ^ 7 + x ^ 5 + x ^ 4 + x ^ 2 + x + 1\]

这和0x104C11DB7的二进制表达式是相等的

CRC积累

有人可能会问: 为什么我不能用硬件加速的CRC32原语并且使用最大的64-bit负载来计算一个更大信息的CRC32? 在数学上是怎么解决的?

如果有人使用原始消息\(m_1\),附加的消息\(m_2\),并且l是\(m_2\)消息的长度,这将是:

\[crc32(m_1||m_2) = x ^ {32} * (m_1||m_2) \pmod p \\\\ = x ^ {32} * (m_1 * 2 ^ l + m_2) \pmod p \\\\ = (x ^ {32 + l} * m_1) + (x ^ {32} * m_2) \pmod p \\\\ = (x ^ l * crc32(m_1)) + (x ^ {32} * m_2) \pmod p\]

用另外的话说，你能使用偏移xor来累计crc值到消息，并且像往常一样的递归使用crc32计算。这说明了CRC32能够像往常一样被递增的累加运算用来计算一个长消息，并且为什么对于硬件实现的CRC32来说拥有累加参数是必要的。

实现模p

我倾向于在x86上面实现CRC32,但我们使用的只是一个无进位的乘法原语。我们没有除法或者取模原语可以使用。别着急！我们能够使用Barrett reduction来实现取模!

Barret reduction来自以下观察:

\[a \pmod p = a - \lfloor sa \rfloor p\]

在该处

\[s = \frac{1}{p}\]

在实践中我们能够近似地使用s地乘法运算为:

\[a \pmod p = a - \lfloor (a * \lfloor \frac{2^{96}}{p} \rfloor) >> 96 \rfloor p\]

译者疑问: 为什么不能直接对除数和被除数使用模2除法呢？应该是\(\lfloor \frac{2 ^ {96}}{p} \rfloor\)可以提前算出来

这能够能够有效地准确算出96位a地准确值。在下面使用地记号中，我们将会使用\(u_i\)来代表\(\lfloor \frac{2^i}{p} \rfloor\)。

Note: 我们选择使用\(u_96\)。有许多的实现使用\(u_64\),很遗憾，改格式只给64-bit地值提供足够地精度。这会导致性能下降，因为这会使用许多非必要的乘法。

译者理解: 这里用96位的目的是因为输入数据a实际上是m(x):64bit + zeros:32bit,所以需要留够96位的精度

实现实战

我们现在对有效实现CRC32的数学基础有很深的理解。但是还是存在一些存在的实现点需要考虑。

Bit-reflected data

译者理解: Bit-reflected指的是高位有效，低位无效的一段bit存储空间，一半的计算指令都是低位有效，高位无效的,总的来说这是一种达到如上效果的一种编码方式

正如你发现的AArch64 CRC32指令使用的ASL伪代码,就像许多的CRC实现,它工作在bit-reflected数据上。幸运的是,不进位的乘法是bit order无感的(没有进位!)。然而，需要考虑的是乘法的的结果存在最小的127位上,就像如下的图片描述(bit position是 bit-reflected):

能够通过将结果左移一位的方式来解决。或者，因为我们的常量乘数总是小于64bits,我们能够提前我们的常量乘数移动1位来获得合适的移位后的结果，如下所述:

类似的，乘一个33位的值能获得一个正确的位对齐。另外一个干净的trick是乘一个65位的最低位为0的数(我们采用这种方法!)。

寻找\(U_i\)

因为\(U_i\)是\(\frac{2^i}{p}\),能够直接通过多项式直接找到:

// Hastily written polynomial division function.
int find_mu(int i) {
    uint256_t dividend = uint256_t{1} << i;
    const uint256_t divisor = 0x104C11DB7;
    const int bits_in_divisor = 33;

    uint256_t result = 0;
    int bit = 255;
    while (bit >= 0) {
        if ((dividend & (uint256_t{1} << bit)) != 0) {
            int shift = bit - bits_in_divisor + 1;
            if (shift >= 0) {
                dividend ^= divisor << shift;
                result ^= uint256_t{1} << shift;
            } else {
                dividend ^= divisor >> -shift;
            }
        }
        bit--;
    }

    printf("%s\n", result.str(16).c_str());
}

记住bit-reflect结果如果你的实现需要它。

容易理解的实现

让我们通过简单的例子来理解实现。记住我们在处理bit-reflected数据,所以谨慎。

using 位反转 = bit-reflected

// This implementation only works for 8-, 16-, 32-bit values
uint32_t crc32(uint32_t value, uint32_t accumulator, int bits) {
    __m128i orig = _mm_set_epi64x(0, (uint64_t)(value ^ accumulator) << (32 - bits));
    __m128i tmp = orig;

    // Multiply by mu_{64}
    tmp = _mm_clmulepi64_si128(tmp, _mm_set_epi64x(0, 0x1F7011641), 0x00);
    // Divide by 2^64 (mask away the unnecessary bits)
    tmp = _mm_and_si128(tmp, _mm_set_epi64x(0, 0xFFFFFFFF));
    // Multiply by p (shifted left by 1 for alignment reasons)
    tmp = _mm_clmulepi64_si128(tmp, _mm_set_epi64x(0, 0x1DB710641), 0x00);
    // Subtract original from result
    tmp = _mm_xor_si128(tmp, orig);

    // Extract the 'lower' (in bit-reflected sense) 32 bits
    return (uint32_t)_mm_extract_epi32(tmp, 1);
}

// For 64-bit values
// Just accumulate using two 32-bit operations
uint32_t crc32(uint64_t value, uint32_t accumulator) {
    accumulator = crc32((uint32_t)value, accumulator, 32);
    accumulator = crc32((uint32_t)(value >> 32), accumulator, 32);
    return accumulator;
}

上述是你能过通过阅读Intel’s white-paper on using PCLMULQDQ to compute CRC能想到的。这本质上是白皮书21页上图12的一个准确实现。

有几个可能提升的点:

1. 常量能够共享一个寄存器
2. 如果你用合适的方法对齐输入,那么没有必要masking
3. 如果选择合适的u值,没必要累计计算64-bit的余数
4. 最后的xor对于32bit和64bit的情况完全没必要

最终实现

uint32_t crc32_32(uint32_t value, uint32_t accumulator) {
    __m128i xmm_const = _mm_set_epi64x(0x00000001DB710641, 0xB4E5B025F7011641);
    __m128i xmm_value = _mm_set_epi64x(0, (value ^ accumulator) << 32);

    xmm_value = _mm_clmulepi64_si128(xmm_value, xmm_const, 0x00);
    xmm_value = _mm_clmulepi64_si128(xmm_value, xmm_const, 0x10);
    return _mm_extract_epi32(xmm_value, 2);
}

uint32_t crc32_64(uint64_t value, uint32_t accumulator) {
    __m128i xmm_const = _mm_set_epi64x(0x00000001DB710641, 0xB4E5B025F7011641);
    __m128i xmm_value = _mm_set_epi64x(0, value ^ accumulator);

    xmm_value = _mm_clmulepi64_si128(xmm_value, xmm_const, 0x00);
    xmm_value = _mm_clmulepi64_si128(xmm_value, xmm_const, 0x10);
    return _mm_extract_epi32(xmm_value, 2);
}

就是这样!通过把32bits移动到右边，我们避免了mask。通过如此设计，我们实现了替代4次乘法的两次乘法的64-bit CRC32。对齐刚好在做96位的移位时不用做任何事情。

我们利用\(u_{92}\)的第0位是0,和输入数据的最低32位由0填充来做\(92-bit \mu 65-bit \rightarrow 156-bit\)的乘法。我们不担心溢出因为我们丢弃答案的最低96位。我们利用这个幸福的巧合64位的输出可以完美的用于下一次乘法。我们然后做一个\(64-bit \mu 33-bit\)的乘法,使用’pextrd’指令截取需要的32bits。

获得的结果就是一个64位CRC32的基于两次乘法的实现。