CUDAの擬似乱数ジェネレーターの概要

作業の詳細により、疑似乱数ジェネレーターを使用してGPUでシミュレーションを行うことがしばしば必要です。 その結果、経験が蓄積され、コミュニティと共有することにしました。

それで、発電機の種類、その長所と短所。

1.線形合同ジェネレータ

周期が2 ^ 32の最も単純で非効率なジェネレーター。 モンテカルロなどの本格的なシミュレーションには絶対に適していませんが、実装が非常に簡単であるため、デバッグに使用できます。

__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. Tauswortheジェネレーターの組み合わせ

このジェネレーターは、LCGよりも長い期間を持っていますが、ほとんど同じくらい簡単です。 周期が2 ^ 113のストリームを生成するには、合計27の命令が必要です。 この期間は、排他的ORを介して複数のジェネレーターを組み合わせることにより実現されます。 以下は、1つのスレッドを生成するコードです。

//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つの符号なし整数を必要としますが、その使用法はLCGと根本的に違いはありません。

3.ハイブリッド発電機

このジェネレーターは、2つの以前のジェネレーターを組み合わせることによって取得されます。 2つのジェネレーターの周期が互いに素である場合、生成されるジェネレーターの周期はそれらの周期の積になります。 このジェネレーターは、3つのTauswortheジェネレーターフローと1つの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つの符号なし整数が必要です。 ヒントをくれたmaratyszczaに感謝します。

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

t = sx^(sx<<11);

sx = sy;
sy = sz;
sz = sw;

sw = (sw^(sw>>19))^(t^(t>>8));

return sw;
}


5.メルセンヌ・ツイスター

このジェネレーターは、現時点で最高のものの1つと考えられており、非常に長い周期（2 ^ 19937）を持っていますが、ストリームから1つの要素を生成するには、複雑な順次アルゴリズムが必要です。 また、ジェネレーターはその状態を保存するために大量のメモリを必要とします。これは、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ジェネレーター

このジェネレーターはGPU専用に開発されたもので、周期は2 ^ 1024です。 ジェネレータは共有メモリを使用するため、ブロック内のすべてのスレッドを協調して実行する必要がありますが、これは常に可能とは限りません。 このジェネレーターの非常に良い説明がここに提供されます。

おわりに

提示された各ジェネレータは、さまざまなケースに適しています。 デバッグにLCGを使用できます。より多くの精度とランダム性が必要な場合は、MersenneTwisterを使用します。 GPU上のすべてのスレッドが一貫した方法で実行される場合-ブロック内に動的な遷移がない場合-WarpStandardは非常に適しています。 それ以外の場合はすべて、HybridGeneratorを使用します。