Поиск списков простых чисел с помощью 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. Поможет ли скаттер-инструкция? Кто-нибудь знает, что это используется на графических процессорах для Решета Эратосфена?

Вот самая быстрая версия Решета Эратосфена с использованием OpenMP, о которой я знаю. Есть идеи, как ускорить это с помощью SSE/AVX? http://create.stephan-brumme.com/eratosthenes/

Вот функция, которую я использую, чтобы определить, является ли число простым. Я работаю сразу с восемью простыми числами (на самом деле это будет 4 за один раз, дважды без AVX2). Я использую векторный класс Агнера Фога. Идея состоит в том, что для больших значений маловероятно, что в последовательности будет восемь простых чисел. Если я найду простое число среди восьми, мне придется перебирать результаты последовательно.

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;
}