起初我想指出模操作可以使用浮點數字來實現:
d % p = d - int(float(d)/float(p))*p.
雖然操作的右側部分的量大然後在左邊這部分是可取的,因爲它可以使用SSE/AVX進行矢量化。
使用SSE4.1的實現:
void Func(const int * a, const int * b, size_t size, int p, int * c)
{
__m128 _k = _mm_set1_ps(1.0f/p);
__m128i _p = _mm_set1_epi32(p);
for (size_t i = 0; i < size; i += 4)
{
__m128i _a = _mm_loadu_si128((__m128i*)(a + i));
__m128i _b = _mm_loadu_si128((__m128i*)(b + i));
__m128i _d = _mm_mul_epi32(_a, _b);
__m128i _e = _mm_cvtps_epi32(_mm_mul_ps(_mm_cvtepi32_ps(_d), _k)); // e = int(float(d)/float(p));
__m128i _c = _mm_sub_epi32(_d, _mm_mul_epi32(_e, _p));
_mm_storeu_si128((__m128i*)(c + i), _c);
}
}
使用AVX的實現:
void Func(const int * a, const int * b, size_t size, int p, int * c)
{
__m256 _k = _mm256_set1_ps(1.0f/p);
__m256i _p = _mm256_set1_epi32(p);
for (size_t i = 0; i < size; i += 8)
{
__m256i _a = _mm256_loadu_si128((__m256i*)(a + i));
__m256i _b = _mm256_loadu_si128((__m256i*)(b + i));
__m256i _d = _mm256_mul_epi32(_a, _b);
__m256i _e = _mm256_cvtps_epi32(_mm256_mul_ps(_mm256_cvtepi32_ps(_d), _k)); // e = int(float(d)/float(p));
__m256i _c = _mm256_sub_epi32(_d, _mm256_mul_epi32(_e, _p));
_mm256_storeu_si128((__m256i*)(c + i), _c);
}
}
這是一個素數。 – Michael
你是否反覆使用同一個'p'? 「尺寸」通常有多大?如果它足夠大或者重複的次數足夠多,使用相同的'p',它可能甚至值得JIT編譯一個使用硬編碼向量移位的循環以及[定點乘法inverse](https://stackoverflow.com/問題/ 41183935 /爲什麼 - 不 - GCC使用乘法按一個怪號碼功能於執行整數-迪維)。或者使用http://libdivide.com/來使用沒有JIT的乘法反轉,但是會有更多的開銷(使用imm8計數的'psrlq'是最好的)。但它可能只有SSE2版本,而不是AVX2或AVX512。 –
'a [i] * b [i]'是否曾經溢出?如果是的話,那會好嗎,還是你想要64位結果mod'p'的結果? – chtz