【问题标题】：AVX 4-bit integersAVX 4 位整数
【发布时间】：2017-10-16 02:52:54
【问题描述】：

我需要执行以下操作：

 w[i] = scale * v[i] + point

比例和点是固定的，而v[] 是一个 4 位整数的向量。

我需要为任意输入向量 v[] 计算 w[]，并且我想使用 AVX 内在函数加速该过程。但是，v[i] 是一个 4 位整数的向量。

问题是如何使用内在函数对 4 位整数执行操作？我可以使用 8 位整数并以这种方式执行操作，但有没有办法执行以下操作：

[a,b] + [c,d] = [a+b,c+d]

[a,b] * [c,d] = [a * b,c * d]

（忽略溢出）

使用 AVX 内部函数，其中 [...,...] 是 8 位整数，a,b,c,d 是 4 位整数？

如果可以，能否举一个简短的例子说明它是如何工作的？

【问题讨论】：

结果溢出怎么办？
您可以使用与SWAR 相同的技术，尽管我怀疑效率是否会比仅解包到 8 位并以每个元素 8 位运行更好。
溢出可以忽略，我们会在 4 位后截断。 SWAR 是一个很好的链接，谢谢！
我仍然不清楚“忽略溢出”是什么意思。如果你添加[1,8]+[1,8]，[3,0] 可以吗，还是你期望[2,0] 或[2,F]？还是 w 应该由 8 位元素组成？
@chtz 像往常一样，是的。如果有几个简单的 8 位乘法，那就太好了..

标签： c++ c vectorization intrinsics avx

【解决方案1】：

只是部分答案（仅添加）和伪代码（应该很容易扩展到 AVX2 内在函数）：

uint8_t a, b;          // input containing two nibbles each

uint8_t c = a + b;     // add with (unwanted) carry between nibbles
uint8_t x = a ^ b ^ c; // bits which are result of a carry
x &= 0x10;             // only bit 4 is of interest
c -= x;                // undo carry of lower to upper nibble

如果已知a 或b 的第4 位未设置（即高半字节的最低位），则可以将其排除在x 的计算之外。

至于乘法：如果 scale 对所有乘积都相同，则您可能会通过一些移位和加/减（在必要时屏蔽溢出位）来摆脱困境。否则，恐怕您需要屏蔽每个 16 位字的 4 位，进行操作，最后将它们放在一起。伪代码（没有AVX 8bit乘法，所以需要用16bit字来操作）：

uint16_t m0=0xf, m1=0xf0, m2=0xf00, m3=0xf000; // masks for each nibble

uint16_t a, b; // input containing 4 nibbles each.

uint16_t p0 = (a*b) & m0; // lowest nibble, does not require masking a,b
uint16_t p1 = ((a>>4) * (b&m1)) & m1;
uint16_t p2 = ((a>>8) * (b&m2)) & m2;
uint16_t p3 = ((a>>12)* (b&m3)) & m3;

uint16_t result = p0 | p1 | p2 | p3;  // join results together

【讨论】：

SWAR 添加的好技巧。这是VPADDB 之后的 4 个额外操作（2 个 VPXOR，1 个 VPAND，1 个 VPSUBB）。第一个异或可以与加法并行运行，因此总延迟 = 4c。显而易见的替代方法是屏蔽事物并将高半字节和低半字节或在一起：VPADB+VPAND 用于低半字节（屏蔽输出），但 2xVPAND + VPADB 用于高半字节（屏蔽输入）。然后VPOR合并。延迟 = 3c, uops = 6 总计，或者比 VPADDB 多 5 个。所以你的方法少了 1 uop，但 dep 链长了 1。所以它非常适合循环遍历两个独立的数组。

【解决方案2】：

4 位加法/乘法可以使用AVX2 完成，特别是如果您想将这些计算应用于更大的向量（比如超过 128 个元素）。但是，如果您只想添加 4 个数字，请使用直接标量代码。

我们在如何处理 4 位整数方面做了大量工作，最近我们开发了一个库来处理它Clover: 4-bit Quantized Linear Algebra Library（重点是量化）。代码也是available at GitHub。

正如您仅提到 4 位整数，我假设您指的是有符号整数（即二进制补码），并据此确定我的答案。请注意，处理无符号实际上要简单得多。

我还假设您希望采用包含n 4 位整数的向量int8_t v[n/2]，并生成具有n/2 4 位整数的int8_t v_sum[n/4]。与下面描述相关的所有代码都是available as a gist。

打包/拆包

显然AVX2 不提供对 4 位整数执行加法/乘法的任何指令，因此，您必须求助于给定的 8 位或 16 位指令。处理 4 位算术的第一步是设计如何将 4 位半字节放入更大的 8 位、16 位或 32 位块中的方法。

为了清楚起见，我们假设您想从一个 32 位块中解压缩给定的半字节，该块将多个 4 位有符号值存储到相应的 32 位整数中（下图）。这可以通过两个位移来完成：

逻辑左移用于移动半字节，使其占据 32 位实体的最高 4 位。
算术右移用于将半字节移动到 32 位实体的最低 4 位。

算术右移有符号扩展，用半字节的符号位填充高位 28 位。产生一个 32 位整数，其值与二进制补码 4 位值相同。

打包（上图左侧）的目标是还原解包操作。两次位移可用于将 32 位整数的最低 4 位放置在 32 位实体中的任何位置。

逻辑左移用于移动半字节，使其占据 32 位实体的最高 4 位。
逻辑右移用于将半字节移动到 32 位实体内的某个位置。

第一个将比半字节低阶的位设置为零，第二个将比半字节高阶的位设置为零。然后可以使用按位或运算在 32 位实体中存储多达 8 个半字节。

如何在实践中应用？

假设您有 64 x 32 位整数值存储在 8 个AVX 寄存器__m256i q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8 中。我们还假设每个值都在 [-8, 7] 范围内。如果要将它们打包到一个 64 x 4 位值的 AVX 寄存器中，可以执行以下操作：

//
// Transpose the 8x8 registers
//
_mm256_transpose8_epi32(q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8);
//
// Shift values left
//
q_1 = _mm256_slli_epi32(q_1, 28);
q_2 = _mm256_slli_epi32(q_2, 28);
q_3 = _mm256_slli_epi32(q_3, 28);
q_4 = _mm256_slli_epi32(q_4, 28);
q_5 = _mm256_slli_epi32(q_5, 28);
q_6 = _mm256_slli_epi32(q_6, 28);
q_7 = _mm256_slli_epi32(q_7, 28);
q_8 = _mm256_slli_epi32(q_8, 28);
//
// Shift values right (zero-extend)
//
q_1 = _mm256_srli_epi32(q_1, 7 * 4);
q_2 = _mm256_srli_epi32(q_2, 6 * 4);
q_3 = _mm256_srli_epi32(q_3, 5 * 4);
q_4 = _mm256_srli_epi32(q_4, 4 * 4);
q_5 = _mm256_srli_epi32(q_5, 3 * 4);
q_6 = _mm256_srli_epi32(q_6, 2 * 4);
q_7 = _mm256_srli_epi32(q_7, 1 * 4);
q_8 = _mm256_srli_epi32(q_8, 0 * 4);
//
// Pack together
//
__m256i t1 = _mm256_or_si256(q_1, q_2);
__m256i t2 = _mm256_or_si256(q_3, q_4);
__m256i t3 = _mm256_or_si256(q_5, q_6);
__m256i t4 = _mm256_or_si256(q_7, q_8);
__m256i t5 = _mm256_or_si256(t1, t2);
__m256i t6 = _mm256_or_si256(t3, t4);
__m256i t7 = _mm256_or_si256(t5, t6);

班次通常需要 1 个周期的吞吐量和 1 个周期的延迟，因此您可以假设这实际上非常便宜。如果您必须处理无符号 4 位值，则可以一起跳过左移。

要反转该过程，您可以应用相同的方法。假设您已将 64 个 4 位值加载到单个 AVX 寄存器 __m256i qu_64 中。为了产生 64 x 32 位整数__m256i q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8，您可以执行以下操作：

//
// Shift values left
//
const __m256i qu_1 = _mm256_slli_epi32(qu_64, 4 * 7);
const __m256i qu_2 = _mm256_slli_epi32(qu_64, 4 * 6);
const __m256i qu_3 = _mm256_slli_epi32(qu_64, 4 * 5);
const __m256i qu_4 = _mm256_slli_epi32(qu_64, 4 * 4);
const __m256i qu_5 = _mm256_slli_epi32(qu_64, 4 * 3);
const __m256i qu_6 = _mm256_slli_epi32(qu_64, 4 * 2);
const __m256i qu_7 = _mm256_slli_epi32(qu_64, 4 * 1);
const __m256i qu_8 = _mm256_slli_epi32(qu_64, 4 * 0);
//
// Shift values right (sign-extent) and obtain 8x8
// 32-bit values
//
__m256i q_1 = _mm256_srai_epi32(qu_1, 28);
__m256i q_2 = _mm256_srai_epi32(qu_2, 28);
__m256i q_3 = _mm256_srai_epi32(qu_3, 28);
__m256i q_4 = _mm256_srai_epi32(qu_4, 28);
__m256i q_5 = _mm256_srai_epi32(qu_5, 28);
__m256i q_6 = _mm256_srai_epi32(qu_6, 28);
__m256i q_7 = _mm256_srai_epi32(qu_7, 28);
__m256i q_8 = _mm256_srai_epi32(qu_8, 28);
//
// Transpose the 8x8 values
//
_mm256_transpose8_epi32(q_1, q_2, q_3, q_4, q_5, q_6, q_7, q_8);

如果处理无符号 4 位，则可以一起跳过右移 (_mm256_srai_epi32)，而不是左移，我们可以执行左逻辑移位 (_mm256_srli_epi32 )。

要查看更多详细信息，请查看gist here。

添加奇数和偶数 4 位条目

假设您使用AVX从向量加载：

const __m256i qv = _mm256_loadu_si256( ... );

现在，我们可以轻松提取奇数和偶数部分。如果AVX2 中有 8 位移位，生活会容易得多，但没有，所以我们必须处理 16 位移位：

const __m256i hi_mask_08   = _mm256_set1_epi8(-16);
const __m256i qv_odd_dirty = _mm256_slli_epi16(qv, 4);
const __m256i qv_odd_shift = _mm256_and_si256(hi_mask_08, qv_odd_dirty);
const __m256i qv_evn_shift = _mm256_and_si256(hi_mask_08, qv);

此时，您在两个AVX 寄存器中基本上分离了奇数和偶数半字节，这些寄存器将它们的值保存在高 4 位（即范围 [-8 * 2^4， 7 * 2^4]）。即使在处理无符号 4 位值时，该过程也是相同的。现在是添加值的时候了。

const __m256i qv_sum_shift = _mm256_add_epi8(qv_odd_shift, qv_evn_shift);

这适用于有符号和无符号，因为二进制加法适用于二进制补码。但是，如果您想避免上溢或下溢，您还可以考虑在 AVX 中已经支持饱和的加法（对于有符号和无符号）：

__m256i _mm256_adds_epi8 (__m256i a, __m256i b)
__m256i _mm256_adds_epu8 (__m256i a, __m256i b)

qv_sum_shift 将在 [-8 * 2^4, 7 * 2^4] 范围内。要将其设置为正确的值，我们需要将其移回（注意如果qv_sum 必须是无符号的，我们可以使用_mm256_srli_epi16 代替）：

const __m256i qv_sum = _mm256_srai_epi16(qv_sum_shift, 4);

求和现已完成。根据您的用例，这也可能是程序的结束，假设您想要生成 8 位的内存块作为结果。但是让我们假设你想解决一个更难的问题。让我们假设输出又是一个 4 位元素的向量，具有与输入相同的内存布局。在这种情况下，我们需要将 8 位块打包成 4 位块。但是，问题在于，我们最终会得到 32 个值（即向量大小的一半），而不是 64 个值。

从这一点来看，有两种选择。我们要么在向量中向前看，处理 128 x 4 位的值，因此我们产生 64 x 4 位的值。或者我们恢复到 SSE，处理 32 x 4 位值。无论哪种方式，将 8 位块打包成 4 位块的最快方法是使用 vpackuswb（或 packuswb 用于 SSE）指令：

__m256i _mm256_packus_epi16 (__m256i a, __m256i b)

该指令使用无符号饱和将压缩的 16 位整数从 a 和 b 转换为压缩的 8 位整数，并将结果存储在 dst 中。这意味着我们必须交错奇数和偶数 4 位值，以便它们驻留在 16 位内存块的 8 个低位中。我们可以进行如下操作：

const __m256i lo_mask_16 = _mm256_set1_epi16(0x0F);
const __m256i hi_mask_16 = _mm256_set1_epi16(0xF0);

const __m256i qv_sum_lo       = _mm256_and_si256(lo_mask_16, qv_sum);
const __m256i qv_sum_hi_dirty = _mm256_srli_epi16(qv_sum_shift, 8);
const __m256i qv_sum_hi       = _mm256_and_si256(hi_mask_16, qv_sum_hi_dirty);

const __m256i qv_sum_16       = _mm256_or_si256(qv_sum_lo, qv_sum_hi);

有符号和无符号 4 位值的过程相同。现在，qv_sum_16 包含两个连续的 4 位值，存储在 16 位内存块的低位中。假设我们从下一次迭代中获得了qv_sum_16（称为qv_sum_16_next），我们可以将所有内容打包为：

const __m256i qv_sum_pack = _mm256_packus_epi16(qv_sum_16, qv_sum_16_next);
const __m256i result      = _mm256_permute4x64_epi64(qv_sum_pack, 0xD8);

或者，如果我们只想生成 32 x 4 位的值，我们可以执行以下操作：

const __m128i lo = _mm256_extractf128_si256(qv_sum_16, 0);
const __m128i hi = _mm256_extractf128_si256(qv_sum_16, 1);
const __m256i result = _mm_packus_epi16(lo, hi)

把它们放在一起

假设有符号的半字节，向量大小n，使得n大于128个元素并且是128的倍数，我们可以执行奇偶加法，产生n/2元素如下：

void add_odd_even(uint64_t n, int8_t * v, int8_t * r)
{
    //
    // Make sure that the vector size that is a multiple of 128
    //
    assert(n % 128 == 0);
    const uint64_t blocks = n / 64;
    //
    // Define constants that will be used for masking operations
    //
    const __m256i hi_mask_08 = _mm256_set1_epi8(-16);
    const __m256i lo_mask_16 = _mm256_set1_epi16(0x0F);
    const __m256i hi_mask_16 = _mm256_set1_epi16(0xF0);

    for (uint64_t b = 0; b < blocks; b += 2) {
        //
        // Calculate the offsets
        //
        const uint64_t offset0 = b * 32;
        const uint64_t offset1 = b * 32 + 32;
        const uint64_t offset2 = b * 32 / 2;
        //
        // Load 128 values in two AVX registers. Each register will
        // contain 64 x 4-bit values in the range [-8, 7].
        //
        const __m256i qv_1 = _mm256_loadu_si256((__m256i *) (v + offset0));
        const __m256i qv_2 = _mm256_loadu_si256((__m256i *) (v + offset1));
        //
        // Extract the odd and the even parts. The values will be split in
        // two registers qv_odd_shift and qv_evn_shift, each of them having
        // 32 x 8-bit values, such that each value is multiplied by 2^4
        // and resides in the range [-8 * 2^4, 7 * 2^4]
        //
        const __m256i qv_odd_dirty_1 = _mm256_slli_epi16(qv_1, 4);
        const __m256i qv_odd_shift_1 = _mm256_and_si256(hi_mask_08, qv_odd_dirty_1);
        const __m256i qv_evn_shift_1 = _mm256_and_si256(hi_mask_08, qv_1);
        const __m256i qv_odd_dirty_2 = _mm256_slli_epi16(qv_2, 4);
        const __m256i qv_odd_shift_2 = _mm256_and_si256(hi_mask_08, qv_odd_dirty_2);
        const __m256i qv_evn_shift_2 = _mm256_and_si256(hi_mask_08, qv_2);
        //
        // Perform addition. In case of overflows / underflows, behaviour
        // is undefined. Values are still in the range [-8 * 2^4, 7 * 2^4].
        //
        const __m256i qv_sum_shift_1 = _mm256_add_epi8(qv_odd_shift_1, qv_evn_shift_1);
        const __m256i qv_sum_shift_2 = _mm256_add_epi8(qv_odd_shift_2, qv_evn_shift_2);
        //
        // Divide by 2^4. At this point in time, each of the two AVX registers holds
        // 32 x 8-bit values that are in the range of [-8, 7]. Summation is complete.
        //
        const __m256i qv_sum_1 = _mm256_srai_epi16(qv_sum_shift_1, 4);
        const __m256i qv_sum_2 = _mm256_srai_epi16(qv_sum_shift_2, 4);
        //
        // Now, we want to take the even numbers of the 32 x 4-bit register, and
        // store them in the high-bits of the odd numbers. We do this with
        // left shifts that extend in zero, and 16-bit masks. This operation
        // results in two registers qv_sum_lo and qv_sum_hi that hold 32
        // values. However, each consecutive 4-bit values reside in the
        // low-bits of a 16-bit chunk.
        //
        const __m256i qv_sum_1_lo       = _mm256_and_si256(lo_mask_16, qv_sum_1);
        const __m256i qv_sum_1_hi_dirty = _mm256_srli_epi16(qv_sum_shift_1, 8);
        const __m256i qv_sum_1_hi       = _mm256_and_si256(hi_mask_16, qv_sum_1_hi_dirty);
        const __m256i qv_sum_2_lo       = _mm256_and_si256(lo_mask_16, qv_sum_2);
        const __m256i qv_sum_2_hi_dirty = _mm256_srli_epi16(qv_sum_shift_2, 8);
        const __m256i qv_sum_2_hi       = _mm256_and_si256(hi_mask_16, qv_sum_2_hi_dirty);
        const __m256i qv_sum_16_1       = _mm256_or_si256(qv_sum_1_lo, qv_sum_1_hi);
        const __m256i qv_sum_16_2       = _mm256_or_si256(qv_sum_2_lo, qv_sum_2_hi);
        //
        // Pack the two registers of 32 x 4-bit values, into a single one having
        // 64 x 4-bit values. Use the unsigned version, to avoid saturation.
        //
        const __m256i qv_sum_pack = _mm256_packus_epi16(qv_sum_16_1, qv_sum_16_2);
        //
        // Interleave the 64-bit chunks.
        //
        const __m256i qv_sum = _mm256_permute4x64_epi64(qv_sum_pack, 0xD8);
        //
        // Store the result
        //
        _mm256_storeu_si256((__m256i *)(r + offset2), qv_sum);
    }
}

此代码的独立测试器和验证器是available in the gist here。

奇数和偶数 4 位条目相乘

对于奇偶条目的乘法，我们可以使用与上述相同的策略将 4 位提取到更大的块中。

AVX2 不提供 8 位乘法，只有 16 位。但是，我们可以按照Agner Fog's C++ vector class library中实现的方法实现8位乘法：

static inline Vec32c operator * (Vec32c const & a, Vec32c const & b) {
    // There is no 8-bit multiply in SSE2. Split into two 16-bit multiplies
    __m256i aodd    = _mm256_srli_epi16(a,8);         // odd numbered elements of a
    __m256i bodd    = _mm256_srli_epi16(b,8);         // odd numbered elements of b
    __m256i muleven = _mm256_mullo_epi16(a,b);        // product of even numbered elements
    __m256i mulodd  = _mm256_mullo_epi16(aodd,bodd);  // product of odd  numbered elements
            mulodd  = _mm256_slli_epi16(mulodd,8);    // put odd numbered elements back in place
    __m256i mask    = _mm256_set1_epi32(0x00FF00FF);  // mask for even positions
    __m256i product = selectb(mask,muleven,mulodd);   // interleave even and odd
    return product;
}

不过，我建议先将半字节提取到 16 位块中，然后使用 _mm256_mullo_epi16 以避免执行不必要的移位。

【讨论】：

来自 cmets 的值可以是 8，所以我认为它是无符号的
如果您仔细阅读说明，您会注意到我还解释了仅支持无符号的简化。

【解决方案3】：

对于w[i]=v[i] * a + b 中的固定a、b，您可以简单地使用查找表w_0_3 = _mm_shuffle_epi8(LUT_03, input) 作为 LSB。将输入拆分为偶数和奇数半字节，奇数 LUT 预移位 4。

auto a = input & 15; // per element
auto b = (input >> 4) & 15; // shift as 16 bits
return LUTA[a] | LUTB[b];

如果有的话，如何动态生成这些 LUT 是另一个问题。

【讨论】：