向量化模算术

Question

我正在尝试编写一些相当快的分量向量加法代码。我使用（有符号，我相信）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）。

我的基准是从 LINE1 中删除 %q。对于 dim=11221184 的100次（迭代）函数调用需要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秒，比没有任何向量化时要慢！这真的是在这里可以做到的最好吗？

- crockeea

---

如果你知道输入已经完全简化，那么最好使用比较和条件减法。 - Mysticial

@PaulR 在运行时，q应该保持（基本）不变，但在编译时是未知的。这有什么优点？ - crockeea

有趣的是，条件减法只花费了1.9秒，这可能是合理的，但ICC 没有向量化。我不知道它为什么这么快。 - crockeea

@Eric，你可以使用SIMD进行条件运算。比较指令会返回一个由全部为0或1的向量组成的结果，你可以将其与另一个值相与并从目标中减去。 - Mysticial

@Eric，根据Mystical的评论，我在我的答案中添加了SSE代码。 - Z boson

Answer 1

既不SSE2，也不AVX2有整数除法指令。英特尔称SVML函数为内部函数是虚伪的，因为其中许多都是复杂的函数，映射到多个指令，而不仅仅是几个。

有一种方法可以使用SSE2或AVX2进行更快的除法（和模除）。请参阅这篇论文 Improved division by invariant integers. 基本上，您需要预先计算一个除数，然后进行乘法运算。预先计算除数需要时间，但对于您代码中某个值的dim，该方法应该是最佳的。我在这里 SSE integer division? 更详细地描述了这种方法。我还在一个质数查找器中成功地实现了这种方法 Finding lists of prime numbers with SIMD - SSE/AVX

Agner Fog在他的Vector Class中使用了该论文中描述的方法实现了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]);
    }   
}

- Z boson