RANDOM, SSE, SSE2

Всё чудесатее и чудесатее. Давеча профайлер показал на rand() в партиклах. В паритклах его много и он там большей частью float rand(), либо vec3 rand(), либо таки vec4 rand(). Дизассемблер немедленно сказал, что стандартный rand() не инлайнится.

В партиклах не нужен “хороший” rand. Он там нужен крайне быстрый. Чем быстрее, тем лучше. Гаденький проканает на ура. Сперва чтиво.

http://en.wikipedia.org/wiki/Linear_congruential_generator

Такой генератор – совершенно банальный вариант, причём некачественный. Если качество начинает влиять, то лучше использовать другой тип генераторов.

http://en.wikipedia.org/wiki/Mersenne_twister

а именно его SFMT подвид

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html

Потому-что он ориентирован на 128-bit SIMD. Как-бы. Т-е SSE2, но не SSE. Потому что целочисленный, хотя перейти к float-ам нетяжело. Правда у него табличка. Вообщем много разных “но”.

Будем изобретать rand для партиклов. Подубовее. Возьмём оригинальный код для вектора, который отнимает “6.42″ попугаев на тестовом примере


__forceinline __m128 mm_rand_std()
{
return _mm_set_ps (
float(rand())/RAND_MAX
, float(rand())/RAND_MAX
, float(rand())/RAND_MAX
, float(rand())/RAND_MAX
);
}

Начнём оптимизации с простого LCG. Итак


unsigned int g_seed = 7; // какой-то seed
unsigned int irand()
{
g_seed = g_seed * 16807;
return (g_seed>>16)&0x7FFF;
}

то что получилось вполне себе инлайнится. Всё то-же самое с float-ами, в дубовом исполнении будет выглядеть вот так


__forceinline float frandi()
{
g_seed = g_seed * 16807;
return float((g_seed>>16) & 0x7fff) * (1.0f/0x7fff);
}

и соответственно вектор


__forceinline __m128 mm_rand_fpui()
{
return _mm_set_ps ( frandi(), frandi(), frandi(), frandi() );
}

Меряем. “0.76″ попугаев. На этом можно было-бы и остановиться, но перфекционизм не позволяет. Хочется быстрее.

Итак фокус – “сделаем” float. Нужные случайные биты сдвинем в мантиссу, а в экспоненту положим биты так, чтобы получилось 1.случайные_биты. Сперва чтиво.

http://en.wikipedia.org/wiki/IEEE_754

Получаем примерно вот такой код


#define MANTISSA_MASK ((1<<23)-1)
#define EXPONENT_MASK 0x3f800000 /* 1.0f */

__forceinline float frand()
{
g_seed = g_seed * 16807;
unsigned int f_result = (g_seed>>8) & MANTISSA_MASK | EXPONENT_MASK;
float result = *((float *)&f_result);
return result - 1.0f;
}

Меряем. “0.84″ попугая. Стало медленнее, что вообщем-то и не удивительно. Но у полученного кода есть замечательное свойство – его можно распараллелить в его общей части, используя только SSE инструкции.


__declspec(align(16)) unsigned int mm_mantissa_mask[4] =
{ MANTISSA_MASK, MANTISSA_MASK, MANTISSA_MASK, MANTISSA_MASK };
__declspec(align(16)) unsigned int mm_exponent_mask[4] =
{ EXPONENT_MASK, EXPONENT_MASK, EXPONENT_MASK, EXPONENT_MASK };

__m128 MM_MANTISSA_MASK = _mm_load_ps ( (const float *) mm_mantissa_mask );
__m128 MM_EXPONENT_MASK = _mm_load_ps ( (const float *) mm_exponent_mask );

unsigned int g_seed0 = 7;
unsigned int g_seed1 = 9;
unsigned int g_seed2 = 13;
unsigned int g_seed3 = 17;

__forceinline __m128 mm_rand_sse()
{
__m128 f_res;
f_res.m128_u32[0]=g_seed0>>8; g_seed0*=16807;
f_res.m128_u32[1]=g_seed1>>8; g_seed1*=16807;
f_res.m128_u32[2]=g_seed2>>8; g_seed2*=16807;
f_res.m128_u32[3]=g_seed3>>8; g_seed3*=16807;
return _mm_sub_ps(_mm_or_ps(_mm_and_ps(f_res,MM_MANTISSA_MASK),
MM_EXPONENT_MASK),MM_EXPONENT_MASK);
}

Меряем. “0.62″ попугая. Стало быстрее, что тоже не удивительно. 20% ускорение, по сравнению с fpui версией – приятный бонус. Приятно, когда результат получается в 10 раз быстрее оригинального кода.

В коде выше есть несколько тонких моментов.

Мы вычитаем MM_EXPONENT_MASK потому, что она равна 1.0f строго – иначе пришлось бы заводить дополнительную константу.

Код получился быстрее во многом потому, что мы разделили seed-ы, и процессор стал выполнять часть умножений параллельно.

Ну и, наконец, SSE2 реализация. Всё совершенно так-же дубово, только умножения распараллеливаем в явном виде.


__m128i MM_GSEED = _mm_set_epi32 ( 7, 9, 13, 17 );
__m128i MM_MUL = _mm_set_epi32 ( 16807, 16807, 16807, 16807 );
__m128i MM_MASK = _mm_set_epi32 ( 0, 0xffffffff, 0, 0xffffffff );

__forceinline __m128 mm_rand_sse2()
{
__m128i seed = MM_GSEED;
__m128i low = _mm_mul_epu32 ( _mm_shuffle_epi32(seed,_MM_SHUFFLE(0,2,1,3)), MM_MUL );
__m128i hig = _mm_mul_epu32 ( _mm_shuffle_epi32(seed,_MM_SHUFFLE(2,0,3,1)), MM_MUL );
low = _mm_and_si128 ( low, MM_MASK );
hig = _mm_and_si128 ( hig, MM_MASK );
low = _mm_shuffle_epi32(low,_MM_SHUFFLE(1,3,0,2));
hig = _mm_shuffle_epi32(hig,_MM_SHUFFLE(0,2,1,3));
MM_GSEED = _mm_or_si128 ( low, hig );
__m128i e_res = _mm_srai_epi32 ( MM_GSEED, 8 );
__m128 f_res = _mm_load_ps ( (const float *) &e_res );
return _mm_sub_ps(_mm_or_ps(_mm_and_ps(f_res,MM_MANTISSA_MASK),
MM_EXPONENT_MASK),MM_EXPONENT_MASK);
}

Меряем. “0.21″ попугая. Хорошо когда есть SSE2. В 3.5 раза быстрее чем fpui, в 30 раз быстрее оригинального кода.

Кстати, код немножко неоптимальный – можно на один SHUFFLE меньше. Или не на один =)

Updated: убран очевидный баг в SSE2 версии. Написан несколько более оптимальный код.

  • digital_cucumber

    У нас тоже вылезал rand() в профайлере, причем тоже именно в партиклах. Насколько я помню, проблема была даже не в том, что он не инлайнится, а в том что имплементация rand() в multithreaded-CRT хранит текущее состояние генератора в thread local storage, и как раз извлечение его оттуда и отжирало львиную долю времени.

  • GeneralGDA

    Аналогичная проблема была. Лечили кэшэм, сгенерированым обычным rand()’ом по старту приложения. Далее выборка из кэша по циклу.

  • look4awhile

    выборка из кеша по циклу – это хорошо, кстати.
    только это ещё один поток данных из памяти. в партиклах будет мешать всему остальному.

  • IronPeter

    А почему не обычное SSE? Если у тебя 16 бит точности – то она отлично выдирается из fp умножения.

  • look4awhile

    потому что нужны нижние биты, а не верхние.
    впрочем напиши. я не осилил

  • Sergei_am

    А вот когда я писал емитер по геометричной форме (примерно 2D окружность), то некачественность rand()-а бъла очень заметна, когда емитилось много партиклов. Очень заметна. И очень заметно бъло как MT его полечил.

  • _asm

    Да. Ранд на спектруме особенно интересен тем что приходится избавляться от умножений.
    Чтобы получилась максимальная длина последоватльности к.г. должен удовлетворять условию – мультипликатор X любое нечетное число и X mod 4 = 1. Слагаемое – любое нечетное. Так что 5 и 1 вполне подходят в условие.
    А чтобы умножвать на 5, вовсе не обязательно задействовать умножение, умножение на 5 можно заменить сдвигом на два бита влево и сложением. Итого:
    [code]
    #include

    typedef unsigned long dword;

    dword SEED = 0x3A092C2A; // any "kickstart"-number

    dword rnd(){
    _asm
    {
    MOV EAX, [SEED]
    SHL EAX, 2
    ADD EAX, [SEED] ; * 5
    INC EAX ; + 1
    MOV [SEED], EAX
    }
    }

    int main(int argc, char* argv) {
    int i;

    for(i=0; i

  • _asm

    Сцуко. :-) раскройте мне тему форматирования кода в wordpress :)

  • http://devlog.gw-labs.net/ dorfe

    _asm, use da http://www.everfall.com/paste/

    а вообще плагин для кода прикрутить надо бы… :-/

  • Pingback: highly professional scums » Blog Archive » RANDOM SSE2()