4

SIMD を使用して素数のリストを見つける方法について誰かアドバイスがあれば知りたいです。特に、SSE/AVX でこれを行う方法に興味があります。

私が注目している 2 つのアルゴリズムは、試行分割とエラトステネスのふるいです。試行分割で SSE を使用する方法を見つけることができました。ベクトル/スカラー「乗算を使用した不変整数による除算」でうまく機能する除算へのより高速な方法を見つけましたhttp://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.1.2556 毎回素数を見つけて結果を形成し、高速除算を行って保存します。次に除算を行うときは、はるかに速くなります。これを行うと、約 3 倍 (4 倍) のスピードアップが得られます。AVX2 を使用すると、さらに高速になる可能性があります。

ただし、試行分割はエラトステネスの篩よりもはるかに遅く、AVX2 でさえまだ行われていないある種のスキャッター命令を除いて、エラトステネスの篩で SIMD を使用する方法は考えられません。分散指導は役に立ちますか?これが GPU でエラトステネスのふるいに使用されていることを知っている人はいますか?

これは、私が知っている OpenMP を使用したエラトステネスのふるいの最速バージョンです。SSE/AVX でこれをスピードアップする方法はありますか? http://create.stephan-brumme.com/eratosthenes/

これは、数値が素数であるかどうかを判断するために使用する関数です。私は一度に 8 つの素数を操作します (実際には、一度に 4 つで、AVX2 なしで 2 回実行されます)。Agner Fog の vectorclass を使用しています。大きな値の場合、8 つの素数が連続して存在する可能性は低いという考えです。8 つの中に素数が見つかった場合は、結果を順番にループする必要があります。

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

ここで、8 つの最有力候補を設定します。これを行うには、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;
}
4

0 に答える 0