CUDA教程三、cuda的随机数生成

上一篇我们介绍了runtime库中的一些函数，接下来我们来介绍cuda随机数的生成。

---

回顾

cuda将函数与变量根据其所在位置，分割成两部分。其中主机端（host）的函数与变量可以互相自由调用，设备端（device）的函数与变量也可自由调用，不过设备端有一种特殊的函数——__kernel__函数（核函数），这是主机端调用设备端函数的唯一方法。

核函数的调用需要运行时参数，host端和device端都是如此。

对于变量的访问，host端可以修改和读取__device__和__constant__，却无法访问__shared__；device端可以读取上述三种，却只能修改__device__和__shared__变量。至于host端变量，自然只有host端函数才能修改和读取，device端无法访问。

问题的提出

神经网络的初始化、游戏NPC的随机生成等，都需要用到随机数。cuda的访问控制，给了随机数生成以一定的问题，host端当然可以用cstdlib中的rand()函数生成随机数，但设备端如何使用这些随机数？每调用一次rand()，就通过cudaMemcpy传递给显存吗？显然不是，这会消耗太多I/O时间；先调用n次，一次性传到device中吗？虽然可行，但并不是最优解。能否用一种方法，让主机仅仅传给设备一个信号，使得多个随机数在device端被生成，这样就能避免过多的I/O带来的耗时问题；或者调用一个设备端函数，使得设备端可以根据需要动态生成随机数。于是curand库和curand_kernel库进入了我们的视野。

curand库——host端的随机数生成

curand库使用了我们的第一种思路，让主机向设备传递信号，指定随机数的生成数量，生成若干随机数以待使用。

随机数生成器

curand提供了多种不同的随机数生成器：

/**
 * CURAND generator types
 */
enum curandRngType {
    CURAND_RNG_TEST = 0,
    CURAND_RNG_PSEUDO_DEFAULT = 100,           ///< Default pseudorandom generator
    CURAND_RNG_PSEUDO_XORWOW = 101,            ///< XORWOW pseudorandom generator
    CURAND_RNG_PSEUDO_MRG32K3A = 121,          ///< MRG32k3a pseudorandom generator
    CURAND_RNG_PSEUDO_MTGP32 = 141,            ///< Mersenne Twister MTGP32 pseudorandom generator
    CURAND_RNG_PSEUDO_MT19937 = 142,           ///< Mersenne Twister MT19937 pseudorandom generator
    CURAND_RNG_PSEUDO_PHILOX4_32_10 = 161,     ///< PHILOX-4x32-10 pseudorandom generator
    CURAND_RNG_QUASI_DEFAULT = 200,            ///< Default quasirandom generator
    CURAND_RNG_QUASI_SOBOL32 = 201,            ///< Sobol32 quasirandom generator
    CURAND_RNG_QUASI_SCRAMBLED_SOBOL32 = 202,  ///< Scrambled Sobol32 quasirandom generator
    CURAND_RNG_QUASI_SOBOL64 = 203,            ///< Sobol64 quasirandom generator
    CURAND_RNG_QUASI_SCRAMBLED_SOBOL64 = 204   ///< Scrambled Sobol64 quasirandom generator
};
typedef enum curandRngType curandRngType_t;

许多英语比我好的同学，通过看注释就能看明白，枚举值大于等于100且小于200的这六个随机数发生器是pseudorandom（伪随机数），而枚举值大于等于200的这五个随机数发生器则是quasirandom（真随机数）。

学过高中数学，我们知道，根据算法得到的随机数为伪随机数，而根据现实中的不依赖算法的随机事件得到的随机数为真随机数。换句话说，伪随机数可以被破解，而真随机数是无法预测的。

比如cstdlib库中的rand函数是这样实现的：

unsigned int _Status = 0;
void __cdecl srand(unsigned int _Seed) {
    _Status = _Seed;
}
int __cdecl rand(void) {
    _Status = 214013u * _Status + 2531011u;
    return (int)(_Status >> 16u & 0x7fffu);
}

很明显这是通过算法模拟出来的随机，当然是伪随机数。

而模拟电视在无信号的时候出现的无序噪点，每一个噪点的值都是真随机数生成的。

curand库生成伪随机数的算法，也是通过算法得到的，如XORWOW算法则是NVIDIA根据XORSHIFT算法改进而成的。而真随机数，则是直接对显卡电流等实际数据的采样，再加上简单的计算处理得到的。在密码学加密中，随机数尽量由真随机数生成器生成，防止敌人预测出序列，轻松破解。

随机数的生成

生成随机数的方法也很简单，首先要初始化随机数生成器，然后调用“传递信号”的函数，让device端生成若干随机数。device端使用时只需知道随机数数组的指针即可。

比如通过下面的操作：

curandGenerator_t gen;        //随机数生成器
curandStatus_t cst;        //curand库函数返回值暂存
cudaError_t cet;        //runtime库函数返回值暂存
unsigned int* arr;        //随机数存储位置
const size_t N = 128;

//初始化生成器
cst = curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_XORWOW);
assert(CURAND_STATUS_SUCCESS == cst);
//动态申请内存
cet = cudaMalloc(&arr, sizeof(unsigned int) * N);
assert(cudaSuccess == cet);
//将N个随机数生成在数组arr中
cst = curandGenerate(gen, arr, N);
assert(CURAND_STATUS_SUCCESS == cst);

本例中使用的是XORWOW生成器，同理我们也可以用其他生成器，包括真随机数生成器，生成其他类型的数、其他类型的分布，同样也可以简略来写：

curandGenerator_t gen;        //随机数生成器
unsigned long long* uint64_arr;    //64位无符号整数数组
double* double_arr;        //64位浮点数数组
const size_t N = 128;

//初始化生成器
assert(CURAND_STATUS_SUCCESS == curandCreateGenerator(&gen, CURAND_RNG_QUASI_SOBOL64));
//动态申请内存
assert(cudaSuccess == cudaMalloc(&uint64_arr, sizeof(unsigned long long) * N));
assert(cudaSuccess == cudaMalloc(&double_arr, sizeof(double) * N));
//将随机数生成在两个数组中
    //生成unsigned long long类型的若干个随机数
assert(CURAND_STATUS_SUCCESS == curandGenerateLongLong(gen, uint64_arr, N));
    //生成符合μ=0、σ²=1标准正态分布的double类型的若干个随机数
assert(CURAND_STATUS_SUCCESS == curandGenerateNormalDouble(gen, double_arr, N, 0.0, 1.0));

除此之外，curand库依然有很多生成不同类型随机数的函数，其中常用的函数原型如下：

curandStatus_t CURANDAPI curandGenerate ( curandGenerator_t generator, unsigned int* outputPtr, size_t num );
//Generate 32-bit pseudo or quasirandom numbers.

curandStatus_t CURANDAPI curandGenerateLogNormal ( curandGenerator_t generator, float* outputPtr, size_t n, float  mean, float  stddev );
//Generate log-normally distributed floats.

curandStatus_t CURANDAPI curandGenerateLogNormalDouble ( curandGenerator_t generator, double* outputPtr, size_t n, double  mean, double  stddev );
//Generate log-normally distributed doubles.

curandStatus_t CURANDAPI curandGenerateLongLong ( curandGenerator_t generator, unsigned long long* outputPtr, size_t num );
//Generate 64-bit quasirandom numbers.

curandStatus_t CURANDAPI curandGenerateNormal ( curandGenerator_t generator, float* outputPtr, size_t n, float  mean, float  stddev );
//Generate normally distributed doubles.

curandStatus_t CURANDAPI curandGenerateNormalDouble ( curandGenerator_t generator, double* outputPtr, size_t n, double  mean, double  stddev );
//Generate normally distributed doubles.

curandStatus_t CURANDAPI curandGeneratePoisson ( curandGenerator_t generator, unsigned int* outputPtr, size_t n, double  lambda );
//Generate Poisson-distributed unsigned ints.

curandStatus_t CURANDAPI curandGenerateUniform ( curandGenerator_t generator, float* outputPtr, size_t num );
//Generate uniformly distributed floats.

curandStatus_t CURANDAPI curandGenerateUniformDouble ( curandGenerator_t generator, double* outputPtr, size_t num );
//Generate uniformly distributed doubles.

编译错误的处理

需要注意的是，用VS直接创建的cuda工程，直接编译上述代码很可能出现一些编译错误：

error LNK2019: 无法解析的外部符号 curandCreateGenerator，函数 main 中引用了该符号
error LNK2019: 无法解析的外部符号 curandGenerate，函数 main 中引用了该符号

有编程经验的朋友们看到LNK就知道这是链接器报错，肯定是少链接了一些文件。没错，只需要在项目-属性-链接器-输入-附加依赖项中添加一项"curand.lib"即可。

curand_kernel库——device端的动态随机数生成

也行有人觉得上述方法太过呆板，能生成多少随机数必须在调用device端函数前指定，实际工程中可能会生成过多的随机数，占用更多内存，也有可能生成的随机数不够，以致访问无效内存或随机数从头循环。

但用第二种思路，使用一些__device__函数动态生成随机数，则会避免这些问题。

可惜的是，国内国外论坛，甚至包括github，很少介绍和使用curand_kernel库，真随机数生成器的初始化方法更是无人介绍，连cuda文档也没有给出例子。

不过也好，我的机会来了——本篇将是历史上第一篇详细介绍curand_kernel库使用方法，并给出使用实例的文章（狗头）。

随机数生成器

curand_kernel提供了八种随机数生成器，并分别封装为结构体：

//伪随机数生成器原型
typedef struct curandStateMtgp32 curandStateMtgp32_t;
typedef struct curandStatePhilox4_32_10 curandStatePhilox4_32_10_t;
typedef struct curandStateXORWOW curandStateXORWOW_t;
typedef struct curandStateMRG32k3a curandStateMRG32k3a_t;
//真随机数生成器原型
typedef struct curandStateSobol32 curandStateSobol32_t;
typedef struct curandStateScrambledSobol32 curandStateScrambledSobol32_t;
typedef struct curandStateSobol64 curandStateSobol64_t;
typedef struct curandStateScrambledSobol64 curandStateScrambledSobol64_t;

名字都很眼熟吧，没错，前四个依旧是伪随机数生成器，后四个依旧是真随机数生成器。除了使用了MTGP32算法的curandStateMtgp32_t生成器外，其他七个生成器都可以用__device__函数——curand_init()初始化。

不过，curand_kernel库的生成器类型并不和curand库的生成器一样是枚举类型，而是定义成了不同的结构体，甚至没使用继承或共用体。这也暗示了不同生成器的初始化需要传入不同的参数。

MRG32K3A、Philox4_32_10与XORWOW生成器的初始化

除去MTGP32生成器外，剩下的三个生成器的初始化方法非常相似：

__device__​ void curand_init ( unsigned long long seed, unsigned long long subsequence, unsigned long long offset, curandStateMRG32k3a_t* state )
//Initialize MRG32k3a state.

__device__​ void curand_init ( unsigned long long seed, unsigned long long subsequence, unsigned long long offset, curandStatePhilox4_32_10_t* state )
//Initialize Philox4_32_10 state.

__device__​ void curand_init ( unsigned long long seed, unsigned long long subsequence, unsigned long long offset, curandStateXORWOW_t* state )
//Initialize XORWOW state.

以XORWOW为例，我们可以在核函数中这么写来初始化多个curandStateXORWOW_t：

__global__ void initialize_generator(curandStateXORWOW_t* states, unsigned long long seed, size_t size) {
    size_t Idx = threadIdx.x + blockDim.x * blockIdx.x;
    if(Idx >= size) return;
    
    curand_init(seed, Idx, 0, &states[Idx]);
}

//调用参数准备
curandStateXORWOW_t* states;
size_t size = 640 * 360;
assert(cudaSuccess == cudaMalloc(&states, sizeof(curandStateXORWOW_t) * size));
//调用运行时参数准备
size_t thread_count = 1024;
size_t block_count = (size - 1) / thread_count + 1;
//调用核函数
initialize_generator<<<block_count, thread_count>>>(states, time(nullptr), size);
cudaDeviceSynchronize();
assert(cudaSuccess == cudaGetLastError());

MTGP32生成器的初始化

而MTGP32生成器的初始化，需要包含curand_mtgp32_host.h库。库中有两个函数：

__host__​ __forceinline__ curandStatus_t curandMakeMTGP32Constants ( const mtgp32_params_fast_t params[], mtgp32_kernel_params_t* p )
//Set up constant parameters for the mtgp32 generator.
__host__ __forceinline__ curandStatus_t CURANDAPI curandMakeMTGP32KernelState ( curandStateMtgp32_t* s, mtgp32_params_fast_t params[], mtgp32_kernel_params_t* k, int  n, unsigned long long seed )
//Set up initial states for the mtgp32 generator.

具体初始化方法如下：

mtgp32_kernel_params* devKernelParams;
curandStateMtgp32_t* states;
size_t size = 128;

assert(cudaSuccess == cudaMalloc(&devKernelParams, sizeof(mtgp32_kernel_params)));
assert(cudaSuccess == cudaMalloc(&states, size * sizeof(curandStateMtgp32_t)));
assert(CURAND_STATUS_SUCCESS == curandMakeMTGP32Constants(mtgp32dc_params_fast_11213, devKernelParams));
assert(CURAND_STATUS_SUCCESS == curandMakeMTGP32KernelState(states, mtgp32dc_params_fast_11213, devKernelParams, size, time(nullptr)));

真随机数生成器的初始化

curand_kernel库中的四个真随机数生成器的初始化都需要用到方向向量：

__device__​ void curand_init ( curandDirectionVectors64_t direction_vectors, unsigned long long scramble_c, unsigned long long offset, curandStateScrambledSobol64_t* state );
//Initialize Scrambled Sobol64 state.
__device__​ void curand_init ( curandDirectionVectors64_t direction_vectors, unsigned long long offset, curandStateSobol64_t* state );
//Initialize Sobol64 state.
__device__​ void curand_init ( curandDirectionVectors32_t direction_vectors, unsigned int  scramble_c, unsigned int  offset, curandStateScrambledSobol32_t* state );
//Initialize Scrambled Sobol32 state.
__device__​ void curand_init ( curandDirectionVectors32_t direction_vectors, unsigned int  offset, curandStateSobol32_t* state );
//Initialize Sobol32 state.

curand库有获取方向向量的函数：

curandStatus_t CURANDAPI curandGetDirectionVectors32 ( curandDirectionVectors32_t* vectors, curandDirectionVectorSet_t set );
//Get direction vectors for 32-bit quasirandom number generation.
curandStatus_t CURANDAPI curandGetDirectionVectors64 ( curandDirectionVectors64_t* vectors, curandDirectionVectorSet_t set );
//Get direction vectors for 64-bit quasirandom number generation.

同时有两个随机数生成器也需要施加干扰。curand库也有获取干扰因子的函数：

curandStatus_t CURANDAPI curandGetScrambleConstants32 ( unsigned int** constants )
//Get scramble constants for 32-bit scrambled Sobol'.
curandStatus_t CURANDAPI curandGetScrambleConstants64 ( unsigned long long** constants )
//Get scramble constants for 64-bit scrambled Sobol'.

具体的生成方法就比较冗长了，这里以Scrambled Sobol 32为例：

__global__ void initialize_generator(curandStateScrambledSobol32_t* states, curandDirectionVectors32_t* devVectors,
                     unsigned int* devScrambleConstants, size_t size) {
    size_t Idx = threadIdx.x + blockDim.x * blockIdx.x;
    if(Idx >= size) return;
    
    curand_init(devVectors[Idx], devScrambleConstants[Idx], Idx * 1024, &states[Idx]);
}

//host端变量
curandDirectionVectors32_t* hostVectors;
unsigned int* hostScrambleConstants;
//device端变量
curandStateScrambledSobol32_t* states;
curandDirectionVectors32_t* devVectors;
unsigned int* devScrambleConstants;

const size_t size = 128;
//获取方向向量和扰动量
assert(CURAND_STATUS_SUCCESS == curandGetDirectionVectors32(&hostVectors, CURAND_SCRAMBLED_DIRECTION_VECTORS_32_JOEKUO6));
assert(CURAND_STATUS_SUCCESS == curandGetScrambleConstants32(&hostScrambleConstants));
//申请device端内存
assert(cudaSuccess == cudaMalloc(&states, size * sizeof(curandStateScrambledSobol32_t)));
assert(cudaSuccess == cudaMalloc(&devVectors, size * sizeof(curandDirectionVectors32_t)));
assert(cudaSuccess == cudaMalloc(&devScrambleConstants, size * sizeof(unsigned int)));
//将获取到的方向向量和扰动量拷贝到device端
assert(cudaSuccess == cudaMemcpy(devVectors, hostVectors, size * sizeof(curandDirectionVectors32_t), cudaMemcpyHostToDevice));
assert(cudaSuccess == cudaMemcpy(devScrambleConstants, hostScrambleConstants, size * sizeof(unsigned int), cudaMemcpyHostToDevice));

size_t thread_count = 256;
size_t block_count = (size - 1) / thread_count + 1;
initialize_generator<<<block_count, thread_count>>>(states, devVectors, devScrambleConstants, size);
cudaDeviceSynchronize();
assert(cudaSuccess == cudaGetLastError());

由于真随机数的初始化没有subsequence参数，所以不同线程随机数的区分要靠offset来完成，因此curand_init的offset参数被我填了Idx * 1024，这个1024不是固定的，只要足够大，使得不同线程生成的随机数互不重合，但又不会太大以至于访问无效内存即可。

生成器的使用

实际使用随机数的时候则直接使用curand函数即可，非常方便。curand函数原型如下：

__device__​ unsigned int     curand ( curandStateMtgp32_t* state )
//Return 32-bits of pseudorandomness from a mtgp32 generator.

__device__​ unsigned long long     curand ( curandStateScrambledSobol64_t* state )
//Return 64-bits of quasirandomness from a scrambled Sobol64 generator.

__device__​ unsigned long long     curand ( curandStateSobol64_t* state )
//Return 64-bits of quasirandomness from a Sobol64 generator.

__device__​ unsigned int     curand ( curandStateScrambledSobol32_t* state )
//Return 32-bits of quasirandomness from a scrambled Sobol32 generator.

__device__​ unsigned int     curand ( curandStateSobol32_t* state )
//Return 32-bits of quasirandomness from a Sobol32 generator.

__device__​ unsigned int     curand ( curandStateMRG32k3a_t* state )
//Return 32-bits of pseudorandomness from an MRG32k3a generator.

__device__​ unsigned int     curand ( curandStatePhilox4_32_10_t* state )
//Return 32-bits of pseudorandomness from an Philox4_32_10 generator.

__device__​ unsigned int     curand ( curandStateXORWOW_t* state )
//Return 32-bits of pseudorandomness from an XORWOW generator.

依旧以XORWOW为例，我们来使用上述代码生成若干可以被2或3整除的随机数：

__global__ void use_generator(curandStateXORWOW_t* states, unsigned int* arr, size_t size) {
    size_t Idx = threadIdx.x + blockDim.x * blockIdx.x;
    if(Idx >= size) return;

    arr[Idx] = curand(&states[Idx]);
    while((arr[Idx] % 2 != 0) && (arr[Idx] % 3 != 0)) {
        arr[Idx] = curand(&states[Idx]);
    }
}

unsigned int* arr;
assert(cudaSuccess == cudaMalloc(&arr, sizeof(unsigned int) * size));
use_generator<<<block_count, thread_count>>>(states, arr, size);
cudaDeviceSynchronize();
assert(cudaSuccess == cudaGetLastError());

除curand()函数外，还有若干其他的生成函数，下例代码我们用MTGP32生成器生成区间内均匀分布的double型变量：

__global__ void use_generator(curandStateMtgp32_t* states, double* arr, size_t size) {
    size_t Idx = threadIdx.x + blockDim.x * blockIdx.x;
    if(Idx >= size) return;

    arr[Idx] = curand_uniform_double(&states[Idx]);
}

double* arr;
assert(cudaSuccess == cudaMalloc(&arr, sizeof(double) * size));
size_t thread_count = 256;
size_t block_count = (size - 1) / thread_count + 1;
use_generator<<<block_count, thread_count>>>(states, arr, size);
cudaDeviceSynchronize();
assert(cudaSuccess == cudaGetLastError());

其他常用的分布函数如下（注意，下文中gen_type可被替换为任何一种curand_kernel库的生成器类型）：

__device__​ float curand_log_normal ( gen_type* state, float  mean, float  stddev );
//Return a log-normally distributed float

__device__​ double curand_log_normal_double ( gen_type* state, double  mean, double  stddev );
//Return a log-normally distributed double

__device__​ float curand_normal ( gen_type* state );
//Return a normally distributed float

__device__​ double curand_normal_double ( gen_type* state );
//Return a normally distributed double

__device__​ unsigned int curand_poisson ( gen_type* state, double lambda );
//Return a Poisson-distributed unsigned int

__device__​ float curand_uniform ( gen_type* state );
//Return a uniformly distributed float

__device__​ double curand_uniform_double ( gen_type* state );
//Return a uniformly distributed double