向量化模运算

我正在尝试编写一些相当快速的分量向量加法代码。我正在使用（我相信是有符号的）64 位整数。

函数是

void addRq (int64_t* a, const int64_t* b, const int32_t dim, const int64_t q) {
    for(int i = 0; i < dim; i++) {
        a[i] = (a[i]+b[i])%q; // LINE1
    }
}

我正在编译icc -std=gnu99 -O3（icc，以便我稍后可以使用 SVML）在 IvyBridge（SSE4.2 和 AVX，但不是 AVX2）上。

我的基线是删除%q来自 LINE1。 100 个（迭代）函数调用dim=11221184需要 1.6 秒。 ICC 自动矢量化 SSE 代码；伟大的。

不过我真的很想做模块化添加。随着%q，ICC 不会自动矢量化代码，并且运行时间为 11.8 秒（！）。即使忽略之前尝试的自动矢量化，这仍然看起来有些过分。

由于我没有 AVX2，因此使用 SSE 进行矢量化需要 SVML，这也许就是 ICC 不自动矢量化的原因。无论如何，这是我对内部循环进行矢量化的尝试：

__m128i qs = _mm_set1_epi64x(q);
for(int i = 0; i < dim; i+=2) {
    __m128i xs = _mm_load_si128((const __m128i*)(a+i));
    __m128i ys = _mm_load_si128((const __m128i*)(b+i));
    __m128i zs = _mm_add_epi64(xs,ys);
    zs = _mm_rem_epi64(zs,qs);
    _mm_store_si128((__m128i*)(a+i),zs);
}

主循环的汇编是：

..B3.4:                         # Preds ..B3.2 ..B3.12
    movdqa    (%r12,%r15,8), %xmm0                          #59.22
    movdqa    %xmm8, %xmm1                                  #60.14
    paddq     (%r14,%r15,8), %xmm0                          #59.22
    call      __svml_i64rem2                                #61.9
    movdqa    %xmm0, (%r12,%r15,8)                          #61.36
    addq      $2, %r15                                      #56.30
    cmpq      %r13, %r15                                    #56.24
    jl        ..B3.4        # Prob 82%                      #56.24

因此代码正在按预期进行矢量化。我知道由于 SVML，我可能无法获得 2 倍的加速，但代码运行时间为 12.5 秒，比根本没有矢量化的情况慢！这真的是这里能做到的最好的事情吗？

SSE2 和 AVX2 都没有整数除法指令。 Intel 将 SVML 函数称为“内在函数”是不诚实的，因为其中许多函数都是复杂的函数，映射到多个指令，而不仅仅是几个指令。

有一种方法可以使用 SSE2 或 AVX2 进行更快的除法（和取模）。请参阅本文改进了不变整数除法。基本上，您预先计算除数，然后进行乘法。预先计算除数需要时间，但对于某个值dim在你的代码中它应该胜出。我在这里更详细地描述了这个方法SSE整数除法？我也在素数查找器中成功实现了这个方法使用 SIMD 查找素数列表 - SSE/AVX

Agner Fog 在他的代码中实现了 32 位（但不是 64 位）除法矢量类使用该论文中描述的方法。如果您想要一些代码，那么这将是一个很好的起点，但您必须将其扩展到 64 位。

编辑：根据 Mysticial 的评论并假设输入已经减少，我为 SSE 制作了一个版本。如果这是在 MSVC 中编译的，那么它需要处于 64 位模式，因为 32 位模式不支持_mm_set1_epi64x。这可以针对 32 位模式进行修复，但我不想这样做。

#ifdef _MSC_VER 
#include <intrin.h>
#endif
#include <nmmintrin.h>                 // SSE4.2
#include <stdint.h>
#include <stdio.h>

void addRq_SSE(int64_t* a, const int64_t* b, const int32_t dim, const int64_t q) {
    __m128i q2 = _mm_set1_epi64x(q);
    __m128i t2 = _mm_sub_epi64(q2,_mm_set1_epi64x(1));
    for(int i = 0; i < dim; i+=2) {
        __m128i a2 = _mm_loadu_si128((__m128i*)&a[i]);
        __m128i b2 = _mm_loadu_si128((__m128i*)&b[i]);
        __m128i c2 = _mm_add_epi64(a2,b2);
        __m128i cmp = _mm_cmpgt_epi64(c2, t2);
        c2 = _mm_sub_epi64(c2, _mm_and_si128(q2,cmp));
        _mm_storeu_si128((__m128i*)&a[i], c2);
    }
}

int main() {
    const int64_t dim = 20;
    int64_t a[dim];
    int64_t b[dim];
    int64_t q = 10;

    for(int i=0; i<dim; i++) {
        a[i] = i%q; b[i] = i%q;
    }
    addRq_SSE(a, b, dim, q);
    for(int i=0; i<dim; i++) {
        printf("%d\n", a[i]);
    }   
}