【发布时间】：2018-12-10 04:47:13
【问题描述】：
我正在编写一个程序来检测素数。一方面是筛选出可能的候选人。我写了一个相当快的程序，但我想我会看看是否有人有更好的想法。我的程序可以使用一些快速收集和分散指令，但我仅限于用于 x86 架构的 AVX2 硬件（我知道 AVX-512 有这些，但我不确定它们有多快）。
#include <stdint.h>
#include <immintrin.h>
#define USE_AVX2
// Sieve the bits in array sieveX for later use
void sieveFactors(uint64_t *sieveX)
{
const uint64_t totalX = 5000000;
#ifdef USE_AVX2
uint64_t indx[4], bits[4];
const __m256i sieveX2 = _mm256_set1_epi64x((uint64_t)(sieveX));
const __m256i total = _mm256_set1_epi64x(totalX - 1);
const __m256i mask = _mm256_set1_epi64x(0x3f);
// Just filling with some typical values (not really constant)
__m256i ans = _mm256_set_epi64x(58, 52, 154, 1);
__m256i ans2 = _mm256_set_epi64x(142, 70, 136, 100);
__m256i sum = _mm256_set_epi64x(201, 213, 219, 237); // 3x primes
__m256i sum2 = _mm256_set_epi64x(201, 213, 219, 237); // This aren't always the same
// Actually algorithm can changes these
__m256i mod1 = _mm256_set1_epi64x(1);
__m256i mod3 = _mm256_set1_epi64x(1);
__m256i mod2, mod4, sum3;
// Sieve until all factors (start under 32-bit threshold) exceed the limit
do {
// Sieve until one of the factors exceeds the limit
do {
// Compiler does a nice job converting these into extracts
*(__m256i *)(&indx[0]) = _mm256_add_epi64(_mm256_srli_epi64(_mm256_andnot_si256(mask, ans), 3), sieveX2);
*(__m256i *)(&bits[0]) = _mm256_sllv_epi64(mod1, _mm256_and_si256(mask, ans));
ans = _mm256_add_epi64(ans, sum);
// Early on these locations can overlap
*(uint64_t *)(indx[0]) |= bits[0];
*(uint64_t *)(indx[1]) |= bits[1];
*(uint64_t *)(indx[2]) |= bits[2];
*(uint64_t *)(indx[3]) |= bits[3];
mod2 = _mm256_sub_epi64(total, ans);
*(__m256i *)(&indx[0]) = _mm256_add_epi64(_mm256_srli_epi64(_mm256_andnot_si256(mask, ans2), 3), sieveX2);
*(__m256i *)(&bits[0]) = _mm256_sllv_epi64(mod3, _mm256_and_si256(mask, ans2));
ans2 = _mm256_add_epi64(ans2, sum2);
// Two types of candidates are being performed at once
*(uint64_t *)(indx[0]) |= bits[0];
*(uint64_t *)(indx[1]) |= bits[1];
*(uint64_t *)(indx[2]) |= bits[2];
*(uint64_t *)(indx[3]) |= bits[3];
mod4 = _mm256_sub_epi64(total, ans2);
} while (!_mm256_movemask_pd(_mm256_castsi256_pd(_mm256_or_si256(mod2, mod4))));
// Remove one factor
mod2 = _mm256_castpd_si256(_mm256_blendv_pd(_mm256_setzero_pd(), _mm256_castsi256_pd(sum), _mm256_castsi256_pd(mod2)));
mod4 = _mm256_castpd_si256(_mm256_blendv_pd(_mm256_setzero_pd(), _mm256_castsi256_pd(sum2), _mm256_castsi256_pd(mod4)));
ans = _mm256_sub_epi64(ans, mod2);
ans2 = _mm256_sub_epi64(ans2, mod4);
sum = _mm256_sub_epi64(sum, mod2);
sum2 = _mm256_sub_epi64(sum2, mod4);
sum3 = _mm256_or_si256(sum, sum2);
} while (!_mm256_testz_si256(sum3, sum3));
#else
// Just some example values (not really constant - compiler will optimize away code incorrectly)
uint64_t cur = 58;
uint64_t cur2 = 142;
uint64_t factor = 67;
if (cur < cur2) {
std::swap(cur, cur2);
}
while (cur < totalX) {
sieveX[cur >> 6] |= (1ULL << (cur & 0x3f));
sieveX[cur2 >> 6] |= (1ULL << (cur2 & 0x3f));
cur += factor;
cur2 += factor;
}
while (cur2 < totalX) {
sieveX[cur2 >> 6] |= (1ULL << (cur2 & 0x3f));
cur2 += factor;
}
#endif
}
请注意，这些位置最初可能会重叠。在循环中短暂停留后，情况并非如此。如果可能的话，我很乐意使用不同的方法。在这部分算法中，大约 82% 的时间都在这个循环中。希望这不会太接近其他已发布的问题。
【问题讨论】：
-
此外，即使使用 AVX512 scatter，可能的重叠（如
cur[0] == cur[1]）意味着您必须检查冲突，以便收集 / SIMD OR / scatter 具有与一次执行一个相同的语义。如果您可以在某个点之后排除这种情况，您可能会使用标量或vpconflictq/ 重试直到那时，然后使用不检查冲突的循环。 -
更新了 gcc 的代码（选项“-mavx2 -O3”），即使我使用的是 MSVC 17。请记住，我“不能”访问 AVX-512 编译器和机器。所以使用 AVX-512 对我没有帮助。我希望能够使用 _mm256_conflict_epi64 和其他高级聚集/分散命令。
-
代码审查：使用
alignas(32) uint64_t *indx[4]，这样您就不必强制转换了。 （并且不要用数组大小 0 声明它，除非这是某种 hack，它可以通过跳过堆栈对齐来获得更有效的编译器输出，因为编译器实际上并没有存储/重新加载。uint64_t indx[0]是错字吗？）我不得不为__m256i sum3;添加一个声明以使其编译（godbolt.org/g/zrSCDa），但是是的，我的回答中的 SIMD 地址计算想法对于 Haswell/Skylake 来说看起来不错。方便（或者不是巧合：P）您已经在向量中拥有正确的掩码。 -
是的，错别字和我删除了所需的变量声明。我只编译了它 - 没有运行它。
-
这太糟糕了标量内存目标
bts不快。bts [rdi], rax将在位串中设置该位，即使它在[rdi]选择的双字之外。 （这种疯狂的 CISC 行为是为什么它并不快！就像 Skylake 上的 10 微秒一样。）
标签： algorithm performance optimization simd avx2