使用 SIMD 查找素数列表 - SSE/AVX

我很好奇是否有人对如何使用 SIMD 查找素数列表有建议。我特别感兴趣如何使用 SSE/AVX 来做到这一点。

我一直在研究的两种算法是试除法和埃拉托斯特尼筛法。 我设法找到一种使用 SSE 进行试除的方法。我找到了一种更快的除法方法，该方法非常适合向量/标量“使用乘法除以不变整数”http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.1.2556每次我找到素数时，我都会形成结果以进行快速除法并保存它们。然后下次我做除法的时候，速度就会快得多。这样做我得到了大约 3 倍的加速（最多 4 倍）。使用 AVX2 也许也会更快。

然而，试除法比埃拉托斯特尼筛法慢得多，我想不出任何方法可以将 SIMD 与埃拉托斯特尼筛法一起使用，除非使用某种分散指令，甚至 AVX2 都还没有。分散指令有帮助吗？有谁知道这个在 GPU 上用于埃拉托斯特尼筛法吗？

这是我所知道的使用 OpenMP 的埃拉托斯特尼筛法的最快版本。有什么想法可以通过 SSE/AVX 来加快速度吗？http://create.stephan-brumme.com/eratosthenes/

这是我用来确定数字是否为素数的函数。我一次对 8 个素数进行操作（实际上，如果没有 AVX2，一次会操作 4 个素数）。我正在使用 Agner Fog 的向量类。这个想法是，对于大值来说，不可能有连续的八个素数。如果我在八个素数中找到一个素数，我必须按顺序循环结果。

inline int is_prime_vec8(Vec8ui num, Divisor_ui_s *buffer, int size) {
    Divisor_ui div = Divisor_ui(buffer[0].m, buffer[0].s1, buffer[0].s2);
    int val = buffer[0].d; 
    Vec8ui cmp = -1;

    for(int i=0; (val*val)<=num[7]; i++) {
        Divisor_ui div = Divisor_ui(buffer[i].m, buffer[i].s1, buffer[i].s2);
        val = buffer[i].d;
        Vec8ui q = num/div; 
        cmp &= (q*val != num);
        int cnt = _mm_movemask_epi8(cmp.get_low()) || _mm_movemask_epi8(cmp.get_high());
        if(cnt == 0) {
            return size;  //0 primes were found
        }
    }
    num &= cmp;  //at least 1 out of 8 values were found to be prime
    int tmp[8];
    num.store(tmp);

    for(int i=0; i<8; i++) {
        if(tmp[i]) {
            set_ui(tmp[i], &buffer[size++]);
        }
    }       
    return size;
}

这是我设置八个主要候选人的地方。我这样做时跳过了 2、3 和 5 的倍数。

int find_primes_vec8(Divisor_ui_s *buffer, const int nmax) {
    int start[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47};
    int size = sizeof(start)/4;

    for(int i=0; i<size; i++) {
        set_ui(start[i], &buffer[i]);
    }

    Vec8ui iv(49, 53, 59, 61, 67, 71, 73, 77);
    size-=3;
    for(int i=49; i<nmax; i+=30) {
        if((i-1)%100000==0) printf("i %d, %f %%\n", i, 100.f*i/(nmax/16));
        size = is_prime_vec8(iv, &buffer[3], size);
        iv += 30;
    }   
    size+=3;

    return size;
}