Pull to refresh

Обзор генераторов псевдослучайных чисел для CUDA

Reading time 5 min
Views 7.4K
По специфике работы часто приходится заниматься симуляциями на GPU с использованием генераторов псевдослучайных чисел. В результате накопился опыт, которым решил и поделиться с сообществом.

Итак, типы генераторов, их достоинства и недостатки.

1. Linear Congruential Generator

Самый простой и неэффективный генератор с периодом всего 2^32. Совершенно не подходит для серьезных симуляций, таких как Monte Carlo, однако можно использовать для отладки, т.к крайне прост в реализации.

__device__ inline unsigned int lcg_rand(unsigned int& seed)
{
seed = seed * 1664525 + 1013904223UL;
return seed;
}


Генератор требует хранить состояние для каждого генератора (потока). Для этого можно использовать массив unsigned int и адресовать его по threadId. Например, так:

__global__ void GenerateRandNums(unsigned int* out, unsigned int* states)
{
const int tid = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int s = states[tid];
out[tid] = lcg_rand(s) + lcg_rand(s);
states[tid] = s;
}


2. Combined Tausworthe Generator

Этот генератор обладает большим периодом чем LCG, однако почти настолько же простой. Он требует всего 27 инструкций, чтобы генерировать поток с периодом 2^113. Такой период достигается путем комбинирования нескольких генераторов посредством exclusive or. Ниже представлен код генерирующий один поток.

//S1,S2,S3,M - константы
__device__ inline unsigned int taus_rand_step(unsigned int& state, int S1, int S2, int S3, int M)
{
unsigned int b = (((state << S1) ^ state) >> S2);
state = (((state & M) << S3) ^ b);
return state;
}


Комбинированный генератор с использованием 4-х независимых потоков.

__device__ unsigned int taus_rand(uint4& state)
{
return
taus_rand_step(state.x, 6, 13, 18, 4294967294UL) ^
taus_rand_step(state.y, 2, 27, 2, 4294967288UL) ^
taus_rand_step(state.z, 13, 21, 7, 4294967280UL) ^
taus_rand_step(state.w, 3, 12, 13, 4294967268UL);

}



Данный генератор требует 4 unsigned int для хранения состояния, однако его использование принципиально не отличается от LCG.

3. Hybrid Generator

Данный генератор получается путем комбинирования 2-х предыдщих генераторов. Если периоды 2-х генераторов взаимно простые, то период результирующего генератора будет произведение их периодов. Данный генератор получается путем комбинирования 3-х потоков Tausworthe генератора и одного потока LCG. Результирующий период будет 2^121.
__device__ unsigned int hybrid_taus(uint4& state)
{
return
taus_rand_step(state.x, 13, 19, 12, 4294967294UL) ^
taus_rand_step(state.y, 2, 25, 4, 4294967288UL) ^
taus_rand_step(state.z, 3, 11, 17, 4294967280UL) ^
lcg_rand(state.w);
}


4. XorShift

Данный генератор обладает периодом 2^128-1 и требует 4 unsigned int для хранения состояния. Спасибо пользователю maratyszcza за совет.

__device__ inline unsigned int rng_xor128(uint4& s)
{
unsigned int t;

t = s.x^(s.x<<11);

s.x = s.y;
s.y = s.z;
s.z = s.w;

s.w = (s.w^(s.w>>19))^(t^(t>>8));

return s.w;
}


5. Mersenne Twister

Данный генератор считается одним из лучших в настоящее время и обладает очень большим периодом — 2^19937, однако для генерации одного элемента из потока требуется выполнение сложного последовательного алгоритма. Также генератор требует много памяти для хранения своего состояния, что может сильно повлиять на производительность GPU, т.к. состояние приходится хранить в глобальной памяти. Однако, в ситуациях когда точность важнее скорости предпостительнее использовать данный генератор. Как и в предыдущих генераторах, у каждого потока есть свой генератор и состояние. Ниже представлен код данного генератора (модифицированный пример из NVIDIA SDK).

#define SEED 999

#define DCMT_SEED 4172
#define MT_RNG_PERIOD 607

typedef struct{
unsigned int matrix_a;
unsigned int mask_b;
unsigned int mask_c;
unsigned int seed;
} mt_struct_stripped;

#define MT_RNG_COUNT (4096*16)
#define MT_MM 9
#define MT_NN 19
#define MT_WMASK 0xFFFFFFFFU
#define MT_UMASK 0xFFFFFFFEU
#define MT_LMASK 0x1U
#define MT_SHIFT0 12
#define MT_SHIFTB 7
#define MT_SHIFTC 15
#define MT_SHIFT1 18

__device__ struct mt_state
{
int iState;
unsigned int mti1;
unsigned int mt[MT_NN];
};

__device__ mt_struct_stripped ds_MT[MT_RNG_COUNT];
static mt_struct_stripped h_MT[MT_RNG_COUNT];
__device__ mt_state ds_mtState[MT_RNG_COUNT];

void InitGPUTwisters(const char* fname, unsigned int seed)
{
FILE *fd = fopen(fname, "rb");
if(!fd)
printf("initMTGPU(): failed to open %s\n", fname);

if( !fread(h_MT, sizeof(h_MT), 1, fd) )
printf("initMTGPU(): failed to load %s\n", fname);

fclose(fd);

mt_struct_stripped *MT = (mt_struct_stripped *)malloc(MT_RNG_COUNT * sizeof(mt_struct_stripped));

for(int i = 0; i < MT_RNG_COUNT; i++)
{
MT[i] = h_MT[i];
seed = lcg_rand(seed);
MT[i].seed = seed;
}
cudaMemcpyToSymbol(ds_MT, MT, sizeof(h_MT));

free(MT);
}

__device__ unsigned int GetRandomIntegerFast(unsigned int tid, mt_state* mts, mt_struct_stripped* config)
{
int iState1, iStateM;
unsigned int mti, mtiM, x;

iState1 = mts->iState + 1;
iStateM = mts->iState + MT_MM;

if(iState1 >= MT_NN)
iState1 -= MT_NN;

if(iStateM >= MT_NN)
iStateM -= MT_NN;

mti = mts->mti1;
mts->mti1 = mts->mt[iState1];
mtiM = mts->mt[iStateM];

x = (mti & MT_UMASK) | (mts->mti1 & MT_LMASK);
x = mtiM ^ (x >> 1) ^ ((x & 1) ? config->matrix_a : 0);
mts->mt[mts->iState] = x;
mts->iState = iState1;

x ^= (x >> MT_SHIFT0);
x ^= (x << MT_SHIFTB) & config->mask_b;
x ^= (x << MT_SHIFTC) & config->mask_c;
x ^= (x >> MT_SHIFT1);

return x;
}

__device__ float GetRandomFloatFast(unsigned int tid, mt_state* mts, mt_struct_stripped* config)
{
return ((float)GetRandomIntegerFast(tid,mts,config) + 1.0f) / 4294967296.0f;
}

__device__ void InitTwisterFast(unsigned int tid, mt_state* mts, mt_struct_stripped* config)
{
mts->mt[0] = config->seed;

for(int i = 1; i < MT_NN; i++)
mts->mt[i] = (1812433253U * (mts->mt[i - 1] ^ (mts->mt[i - 1] >> 30)) + i) & MT_WMASK;

mts->iState = 0;
mts->mti1 = mts->mt[0];
}



Пример использования данного генератора:

__global__ void GenerateUniformNumbers(float* output)
{
const uint tid = gridDim.x*gridDim.y*threadIdx.x + blockIdx.y*gridDim.x + blockIdx.x;

mt_state mts = ds_mtState[tid];
mt_struct_stripped config = ds_MT[tid];

InitTwisterFast(tid, &mts, &config);

output[tid] = GetRandomFloatFast(tid, &mts, &config)

ds_mtState[tid] = mts;
}


6. WarpStandard generator

Данный генератор был разработан специально для GPU и обладает периодом 2^1024. Генератор использует shared memory, таком образом все потоки в блоке должны выполняться согласованно, что не всегда возможно. Очень хорошее описание данного генератора представлено по ссылке

Заключение

Каждый из представленных генераторов подходит для разных случаев. Для отладки можно использовать LCG, если нужна обчень большая точность и randomness я бы использовал MersenneTwister. Если все потоки на GPU выполняются согласовано — нет динамических переходов внутри блоков — очень хорошо подходит WarpStandard. Во всех остальных случаях я использую HybridGenerator.
Tags:
Hubs:
+37
Comments 13
Comments Comments 13

Articles