AVX/SSE версия xorshift128 +

Я стараюсь сделать максимально возможный высококачественный RNG. Прочитав http://xorshift.di.unimi.it/, xorshift128+ кажется хорошим вариантом. Код C

#include <stdint.h>
uint64_t s[ 2 ];

uint64_t next(void) { 
    uint64_t s1 = s[ 0 ];
    const uint64_t s0 = s[ 1 ];
    s[ 0 ] = s0;
    s1 ^= s1 << 23; // a
    return ( s[ 1 ] = ( s1 ^ s0 ^ ( s1 >> 17 ) ^ ( s0 >> 26 ) ) ) + s0; // b, c
}

Я не эксперт SSE/AVX, но мой процессор поддерживает инструкции SSE4.1/SSE4.2/AVX/F16C/FMA3/XOP. Как вы могли использовать их для ускорения этого кода (если вы хотите сделать миллиарды таких случайных чисел) и каков ожидаемый предел этого ускорения на практике?

Ответ 1

XorShift действительно хороший выбор. Это так хорошо, так быстро и требует такого небольшого состояния, что я удивлен, увидев столь малое принятие. Он должен быть стандартным генератором на всех платформах. Я сам реализовал его 8 лет назад, и даже тогда он может генерировать 800 Мбайт/с случайных байтов.

Вы не можете использовать векторные инструкции для ускорения генерации одного случайного числа. В этих немногих инструкциях слишком мало параллелизма на уровне инструкций.

Но вы можете легко ускорить генерацию N чисел, где N - векторный размер вашего набора заданных инструкций. Просто запустите N генераторов параллельно. Сохранять состояние для N генераторов и генерировать N чисел одновременно.

Если клиентский код требует номера по одному, вы можете сохранить буфер из N (или более) номеров. Если буфер пуст, вы заполняете его с помощью векторных инструкций. Если буфер не пуст, вы просто возвращаете следующий номер.

Ответ 2

Для тех, кто может решить этот вопрос, я думаю, что этот код C++ правильно реализует 4 генератора xorshift128plus, работающих параллельно, используя AVX2:

__m256i xorshift128plus_avx2(__m256i &state0, __m256i &state1)
{
    __m256i s1 = state0;
    const __m256i s0 = state1;
    state0 = s0;
    s1 = _mm256_xor_si256(s1, _mm256_slli_epi64(s1, 23));
    state1 = _mm256_xor_si256(_mm256_xor_si256(_mm256_xor_si256(s1, s0),
                                               _mm256_srli_epi64(s1, 18)),
                              _mm256_srli_epi64(s0, 5));
    return _mm256_add_epi64(state1, s0);
}

Скалярная реализация, которую я использовал:

u64 xorshift128plus(u64 &state0, u64 &state1)
{
    u64 s1 = state0;
    const u64 s0 = state1;
    state0 = s0;
    s1 ^= s1 << 23;                              // a
    state1 = s1 ^ s0 ^ (s1 >> 18) ^ (s0 >> 5); // b, c
    return state1 + s0;
}

Это то же самое в бумаге xorshiftplus. Заметим, что константы правого сдвига от исходного вопроса не соответствуют тем, которые указаны в статье.