模乘向量化
Vectorization of modulo multiplication
我有一个功能:
void Func(const int * a, const int * b, size_t size, int p, int * c)
{
for (size_t i = 0; i < size; ++i)
c[i] = (a[i]*b[i])%p;
}
该函数对整数数组进行多次模乘。
所有整数都是正数。
我需要提高它的性能。
我想到了 SSE 和 AVX。但是他们没有向量化模乘法的操作。
或许我错了?
也许有人知道解决这个问题的可能性吗?
首先我要注意的是,模运算可以通过使用浮点数来实现:
d % p = d - int(float(d)/float(p))*p.
尽管右侧部分的运算量比左侧部分大,但此部分更可取,因为它可以使用 SSE/AVX.
进行矢量化
32x32 => 32-bit integer multiplication 的 SSE4.1 实现。请注意,从 FP 转换回整数是通过舍入到最近完成的;如果你想要像 C float->integer conversions 这样的语义,使用向零截断(cvttps_epi32
)。
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_mullo_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_mullo_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_mullo_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_mullo_epi32(_e, _p));
_mm256_storeu_si128((__m256i*)(c + i), _c);
}
}
实际上有一个内在函数正在执行此操作:
_mm256_irem_epi32
https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=_mm256_irem_epi32
我有一个功能:
void Func(const int * a, const int * b, size_t size, int p, int * c)
{
for (size_t i = 0; i < size; ++i)
c[i] = (a[i]*b[i])%p;
}
该函数对整数数组进行多次模乘。 所有整数都是正数。 我需要提高它的性能。
我想到了 SSE 和 AVX。但是他们没有向量化模乘法的操作。 或许我错了?
也许有人知道解决这个问题的可能性吗?
首先我要注意的是,模运算可以通过使用浮点数来实现:
d % p = d - int(float(d)/float(p))*p.
尽管右侧部分的运算量比左侧部分大,但此部分更可取,因为它可以使用 SSE/AVX.
进行矢量化32x32 => 32-bit integer multiplication 的 SSE4.1 实现。请注意,从 FP 转换回整数是通过舍入到最近完成的;如果你想要像 C float->integer conversions 这样的语义,使用向零截断(cvttps_epi32
)。
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_mullo_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_mullo_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_mullo_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_mullo_epi32(_e, _p));
_mm256_storeu_si128((__m256i*)(c + i), _c);
}
}
实际上有一个内在函数正在执行此操作: _mm256_irem_epi32
https://software.intel.com/sites/landingpage/IntrinsicsGuide/#text=_mm256_irem_epi32