逐次で疑似乱数を生成する時、始点からドミノ倒し的に乱数を生成するので、並列化出来ない問題がある。
CUDAはこの問題を解決するAPI[1]を用意しており、<curand.h>、<curand_kernel.h>をincludeすることで使用できる。
CUDAで乱数生成の並列化を行う際は、次のような段取りを取る。
curandState の初期化、シード生成curandState から乱数を生成curandState の更新(2, 3は同時に行われる)curandStateとは乱数を生成するための無秩序な状態ベクトルで、ここから何らかの変換を行うことで様々な乱数を生成する。
CUDAの乱数はcurandState型の状態ベクトルから生成される。
なので、まずは curandState *state の初期化を行わなければならない。
それには
__device__ void curand_init(unsigned long long seed,unsigned long long sequence, unsigned long long offset, curandState *state)
を用いる[1]。
seed はシード値、sequence にはベクトルの要素、*state には状態ベクトルの要素のアドレスを入れる。offset は乱数列をどれだけスキップするかを指定する(一般的には0で良い)。
以上を踏まえて初期化のカーネルを組むと以下のようになる。
#include <cuda.h>
#include <curand_kernel.h>
__global__ void setCurand(unsigned long long seed, curandState *state){
int i_global = threadIdx.x + blockIdx.x * blockDim.x;
curand_init(seed, i_global, 0, &state[i_global]);
}
次に、初期化した*stateを用いて作ることができる基本的な乱数を紹介する。 次に紹介するもの以外の様々な乱数を生成することも出来るがそれはドキュメントを参考にしてほしい[1]。
__device__ unsigned int curand(curandState *state) はunsigned int32型の一様乱数を返す。
__global__ genrand_uint32(unsigned int *result, curandState *state){
uint i_global = threadIdx.x + blockIdx.x*blockDim.x;
//resultに乱数が返り、state[]は更新される。
result[i_global] = curand(&state[i_global]);
}
__device__ float curand_uniform(curandState *state) は、範囲0から1のfloat型の一様乱数を返す。
__global__ genrand_uniform(float *result, curandState *state){
uint i_global = threadIdx.x + blockIdx.x*blockDim.x;
//resultに乱数が返り、state[]は更新される。
result[i_global] = curand_uniform(&state[i_global]);
}
__device__ float curand_normal(curandState * state) は、平均0、分散1の正規乱数を返す。
__global__ genrand_normal(float *result, curandState *state){
uint i_global = threadIdx.x + blockIdx.x*blockDim.x;
//resultに乱数が返り、state[]は更新される。
result[i_global] = curand_normal(&state[i_global]);
}
実際に正規乱数を生成し、平均値と標準偏差を調べてみた。 なお、シードにはメルセンヌ・ツイスタを用いた。
#include <cuda.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <curand.h>
#include <curand_kernel.h>
#include "MT.h"
//スレッド数
#define NT 128
//ブロック数
#define NB 4
//curandStateの初期化
__global__ void setCurand(unsigned long long seed, curandState *state){
uint i_global = threadIdx.x + blockIdx.x*blockDim.x;
curand_init(seed, i_global, 0, &state[i_global]);
}
//一様乱数を返す
__global__ void genrand_kernel(float *result, curandState *state){
uint i_global = threadIdx.x + blockIdx.x*blockDim.x;
result[i_global] = curand_normal(&state[i_global]);
}
int main(){
//Mersenne Twisterの初期化
init_genrand((unsigned long)time(NULL));
//乱数の配列
float *rnd, *rnd_dev;
//av: 平均値, sigma: 標準偏差
double av = 0.;
double sigma = 0.;
//乱数の状態
curandState *state;
rnd = (float*)malloc(NB* NT * sizeof(float));
cudaMalloc((void**)&rnd_dev, NB * NT * sizeof(float));
cudaMalloc((void**)&state, NB * NT * sizeof(curandState));
//Mersenne Twisterをシードに取ってcurandStateを初期化
setCurand<<<NB,NT>>>(genrand_int32(void)), state);
//正規乱数の生成
genrand_kernel<<<NB,NT>>>(rnd_dev,state);
cudaMemcpy(rnd, rnd_dev, NB * NT * sizeof(uint),cudaMemcpyDeviceToHost);
for(uint i = 0; i < NB * NT; i++){
printf("%d, %f\n",i, rnd[i]);
av += rnd[i]/(NB*NT);
}
for(uint i = 0; i < NB * NT; i++){
sigma += (rnd[i]-av)*(rnd[i]-av)/(NB*NT);
}
sigma = sqrt(sigma);
//平均と標準偏差を表示
printf("\nav = %f\nsig = %f\n",av,sigma);
free(rnd);
cudaFree(rnd_dev);
cudaFree(state);
return 0;
}
計算機では熱雑音等を用いない限り、完全な乱数を生成することはできず、周期性がある。
今回の擬似乱数の生成にはXORWOWと呼ばれる手法が用いられている[1]。 状態空間は160bitで周期性は10^48ほどである。 これは一般的なシミュレーションでは使い果たすことのないほどの周期性であり、計算機科学上、乱数として用いて問題ないとされる[2]。