在 GCC 10.3.0 中找不到 _mm256_rem_epu64 内在函数

Question

我尝试用AVX-512 指令重写以下uint64_t 2x2 矩阵乘法，但GCC 10.3 找不到_mm256_rem_epu64 内在函数。

#include <cstdint>
#include <immintrin.h>

constexpr uint32_t LAST_9_DIGITS_DIVIDER = 1000000000;

void multiply(uint64_t f[2][2], uint64_t m[2][2])
{
  uint64_t x = (f[0][0] * m[0][0] + f[0][1] * m[1][0]) % LAST_9_DIGITS_DIVIDER;
  uint64_t y = (f[0][0] * m[0][1] + f[0][1] * m[1][1]) % LAST_9_DIGITS_DIVIDER;
  uint64_t z = (f[1][0] * m[0][0] + f[1][1] * m[1][0]) % LAST_9_DIGITS_DIVIDER;
  uint64_t w = (f[1][0] * m[0][1] + f[1][1] * m[1][1]) % LAST_9_DIGITS_DIVIDER;

  f[0][0] = x;
  f[0][1] = y;
  f[1][0] = z;
  f[1][1] = w;
}

void multiply_simd(uint64_t f[2][2], uint64_t m[2][2])
{
  __m256i v1 = _mm256_set_epi64x(f[0][0], f[0][0], f[1][0], f[1][0]);
  __m256i v2 = _mm256_set_epi64x(m[0][0], m[0][1], m[0][0], m[0][1]);
  __m256i v3 = _mm256_mullo_epi64(v1, v2);

  __m256i v4 = _mm256_set_epi64x(f[0][1], f[0][1], f[1][1], f[1][1]);
  __m256i v5 = _mm256_set_epi64x(m[1][0], m[1][1], m[1][0], m[1][1]);
  __m256i v6 = _mm256_mullo_epi64(v4, v5);

  __m256i v7 = _mm256_add_epi64(v3, v6);
  __m256i div = _mm256_set1_epi64x(LAST_9_DIGITS_DIVIDER);
  __m256i v8 = _mm256_rem_epu64(v7, div);
  _mm256_store_epi64(f, v8);
}

是否有可能以某种方式启用_mm256_rem_epu64，或者如果没有，是否可以通过其他方式使用 SIMD 指令计算提醒？

Answer 1

正如 Peter Cordes 在 cmets 中提到的，_mm256_rem_epu64 是一个 SVML 函数。大多数编译器不支持 SVML； AFAIK 真的只有 ICC 可以，但 clang 也可以配置为使用它。

我知道的唯一其他 SVML 实现是在我的一个项目中，SIMDe。在这种情况下，由于您使用的是 GCC 10.3，the implementation of _mm256_rem_epu64 将使用 vector extensions，因此 the code from SIMDe 将与以下内容基本相同：

#include <immintrin.h>
#include <stdint.h>

typedef uint64_t u64x4 __attribute__((__vector_size__(32)));

__m256i
foo_mm256_rem_epu64(__m256i a, __m256i b) {
    return (__m256i) (((u64x4) a) % ((u64x4) b));
}

在这种情况下，GCC 和 clang 都会对操作进行标量化（参见 Compiler Explorer），因此性能会很差，尤其是考虑到 how slow div 指令。

也就是说，由于您使用的是编译时常量，编译器应该能够用乘法和移位替换除法，因此性能会更好，但我们可以通过使用 @987654327 挤出更多@。

Libdivide 通常在运行时计算魔法值，但libdivide_u64_t 结构非常简单，我们可以跳过libdivide_u64_gen 步骤并在编译时提供结构：

__m256i div_by_1000000000(__m256i a) {
    static const struct libdivide_u64_t d = {
        UINT64_C(1360296554856532783),
        UINT8_C(93)
    };
    return libdivide_u64_do_vec256(a, &d);
}

现在，如果您可以使用 AVX-512VL + AVX-512DQ，则有一个 64 位乘法函数 (_mm256_mullo_epi64)。如果您可以使用它，那可能是正确的方法：

__m256i rem_1000000000(__m256i a) {
    static const struct libdivide_u64_t d = {
        UINT64_C(1360296554856532783),
        UINT8_C(93)
    };
    return
        _mm256_sub_epi64(
            a,
            _mm256_mullo_epi64(
                libdivide_u64_do_vec256(a, &d),
                _mm256_set1_epi64x(1000000000)
            )
        );
}

（或在Compiler Explorer，使用 LLVM-MCA）

如果您没有 AVX-512DQ+VL​​，您可能希望再次使用矢量扩展：

typedef uint64_t u64x4 __attribute__((__vector_size__(32)));

__m256i rem_1000000000(__m256i a) {
    static const struct libdivide_u64_t d = {
        UINT64_C(1360296554856532783),
        UINT8_C(93)
    };
    u64x4 one_billion = { 1000000000, 1000000000, 1000000000, 1000000000 };
    return (__m256i) (
        (
            (u64x4) a) -
            (((u64x4) libdivide_u64_do_vec256(a, &d)) * one_billion
        )
    );
}

(on Compiler Explorer)

所有这些都未经测试，但假设我没有犯任何愚蠢的错误，它应该相对活泼。

如果您真的想摆脱 libdivide 依赖项，您可以自己执行这些操作，但我真的没有理由不使用 libdivide，所以我将把它留给其他人练习。