CUDAで乱数を使いたい!
今回の記事はCUDAで疑似乱数をデバイス上に並列に呼び出したい人に向けたものである。日本語の文献が少ないので、必要な文だけ(おそらく読者にとっても)、ここにまとめることにする。
私の専門はガラスの解析なのだが、最近、GPUを使って分子動力学法をすることにチャレンジしている。Langevin方程式を元にブラウン動力学法を行う際、熱として乱数を用いる必要があるが、逐次で呼び出すわけには行かないので、ちゃんと並列に乱数を呼び出したい。
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でいい)。
以上を踏まえて初期化のカーネルを組むと以下のようになる。
次に、初期化した*stateを用いて作ることができる基本的な乱数を紹介する。
一様乱数(非負整数型)
__device__ uint curand(curandState *state)は32bitのuint型の一様乱数を返す。
一様乱数(0,1]
__device__ float curand_uniform(curandState *state)は、範囲0から1のfloat型の一様乱数を返す。
正規乱数
__device__ float curand_normal(curandState * state)は、平均0、分散1の正規乱数を返す。
実際に使用してみる
実際に正規乱数を生成し、平均値と標準偏差を調べてみよう。なお、シードにはメルセンヌ・ツイスタを用いた。[2, 3]
周期性
PCでは熱雑音等を用いない限り、完全な乱数を生成することはできず、周期性がある。
今回の擬似乱数の生成にはXORWOWと呼ばれる手法が用いられている[1]。状態空間は160bitで周期性は10^48ほどである。これは一般的なシミュレーションでは使い果たすことのないほどの周期性であり、計算機科学上、乱数として用いて問題ないとされる[4]。
References
[1] https://docs.nvidia.com/cuda/curand/device-api-overview.html#device-api-overview
[2] https://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html
[3] https://www.sat.t.u-tokyo.ac.jp/~omi/random_variables_generation.html