1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
#ifndef __RANDOM_H__
#define __RANDOM_H__
#include <curand_kernel.h>
#include "physical_constants.h"
__device__ float
uniform(curandState *s, const float &low, const float &high)
{
return low + curand_uniform(s)*(high-low);
}
__device__ float3
uniform_sphere(curandState *s)
{
float theta = uniform(s, 0.0f, 2*PI);
float u = uniform(s, -1.0f, 1.0f);
float c = sqrtf(1.0f-u*u);
return make_float3(c*cosf(theta), c*sinf(theta), u);
}
// Draw a random sample given a cumulative distribution function
// Assumptions: ncdf >= 2, cdf_y[0] is 0.0, and cdf_y[ncdf-1] is 1.0
__device__ float
sample_cdf(curandState *rng, int ncdf, float *cdf_x, float *cdf_y)
{
float u = curand_uniform(rng);
// Find u in cdf_y with binary search
// list must contain at least 2 elements: 0.0 and 1.0
int lower=0;
int upper=ncdf-1;
while(lower < upper-1) {
int half = (lower+upper) / 2;
if (u < cdf_y[half])
upper = half;
else
lower = half;
}
float frac = (u - cdf_y[lower]) / (cdf_y[upper] - cdf_y[lower]);
return cdf_x[lower] * frac + cdf_x[upper] * (1.0f - frac);
}
extern "C"
{
__global__ void
init_rng(int nthreads, curandState *s, unsigned long long seed,
unsigned long long offset)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
if (id >= nthreads)
return;
curand_init(seed, id, offset, &s[id]);
}
__global__ void
fill_uniform(int nthreads, curandState *s, float *a, float low, float high)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
if (id >= nthreads)
return;
a[id] = uniform(&s[id], low, high);
}
__global__ void
fill_uniform_sphere(int nthreads, curandState *s, float3 *a)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
if (id >= nthreads)
return;
a[id] = uniform_sphere(&s[id]);
}
} // extern "c"
#endif
|
