CUDAでの乱数

CUDAでの乱数†

C for CUDAではC言語のrand関数のような乱数を生成する関数は用意されていません． C言語のrand関数は乱数の質としては問題があるが気軽に使えるので便利である． CUDA SDKにはメルセンヌツイスターのサンプルがあるのでこれを使うこともできるが， ここでは，より簡単に実装できる線形合同法とXORシフト法をデバイス関数として実装してみた．

線形合同法†

C言語のrand関数などで使われている方法．

__device__
unsigned int Rand(unsigned int randx)
{
    randx = randx*1103515245+12345;
    return randx&2147483647;
}

ただし，スレッドごとに異なるrandxを用意して，返値のrandxで更新していく必要がある． また，一度のカーネル関数呼び出し内でスレッドごとに乱数を用いたい場合は， 事前に乱数テーブルをデバイスメモリ(テクスチャメモリやコンスタントメモリがよいかも)に格納しておいて用いる．

XORシフト法†

"Xorshift RNGs"という論文 *1 で提案されている方法． 高速で質もよく，簡単に実装できるのでお勧め．

__device__ static unsigned long xors_x = 123456789;
__device__ static unsigned long xors_y = 362436069;
__device__ static unsigned long xors_z = 521288629;
__device__ static unsigned long xors_w = 88675123;

__device__
unsigned long Xorshift128()
{ 
    unsigned long t; 
    t = (xors_x^(xors_x<<11));
	xors_x = xors_y; xors_y = xors_z; xors_z = xors_w; 
	return ( xors_w = (xors_w^(xors_w>>19))^(t^(t>>8)) ); 
}
__device__
long Xorshift128(long l, long h)
{ 
    unsigned long t; 
    t = (xors_x^(xors_x<<11));
	xors_x = xors_y; xors_y = xors_z; xors_z = xors_w; 
	xors_w = (xors_w^(xors_w>>19))^(t^(t>>8));
	return l+(xors_w%(h-l));
}

__device__
float XorFrand(float l, float h)
{
	return l+(h-l)*(Xorshift128(0, 1000000)/1000000.0f);
}