Fix #35665: cycles CUDA crash after recent changes. This works around a compiler
[blender.git] / intern / cycles / kernel / kernel_random.h
1 /*
2  * Copyright 2011, Blender Foundation.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18
19 #include "kernel_jitter.h"
20
21 CCL_NAMESPACE_BEGIN
22
23 typedef uint RNG;
24
25 #ifdef __SOBOL__
26
27 /* skip initial numbers that are not as well distributed, especially the
28  * first sequence is just 0 everywhere, which can be problematic for e.g.
29  * path termination */
30 #define SOBOL_SKIP 64
31
32 /* High Dimensional Sobol */
33
34 /* van der corput radical inverse */
35 __device uint van_der_corput(uint bits)
36 {
37         bits = (bits << 16) | (bits >> 16);
38         bits = ((bits & 0x00ff00ff) << 8) | ((bits & 0xff00ff00) >> 8);
39         bits = ((bits & 0x0f0f0f0f) << 4) | ((bits & 0xf0f0f0f0) >> 4);
40         bits = ((bits & 0x33333333) << 2) | ((bits & 0xcccccccc) >> 2);
41         bits = ((bits & 0x55555555) << 1) | ((bits & 0xaaaaaaaa) >> 1);
42         return bits;
43 }
44
45 /* sobol radical inverse */
46 __device uint sobol(uint i)
47 {
48         uint r = 0;
49
50         for(uint v = 1U << 31; i; i >>= 1, v ^= v >> 1)
51                 if(i & 1)
52                         r ^= v;
53
54         return r;
55 }
56
57 /* inverse of sobol radical inverse */
58 __device uint sobol_inverse(uint i)
59 {
60         const uint msb = 1U << 31;
61         uint r = 0;
62
63         for(uint v = 1; i; i <<= 1, v ^= v << 1)
64                 if(i & msb)
65                         r ^= v;
66
67         return r;
68 }
69
70 /* multidimensional sobol with generator matrices
71  * dimension 0 and 1 are equal to van_der_corput() and sobol() respectively */
72 __device uint sobol_dimension(KernelGlobals *kg, int index, int dimension)
73 {
74         uint result = 0;
75         uint i = index;
76
77         for(uint j = 0; i; i >>= 1, j++)
78                 if(i & 1)
79                         result ^= kernel_tex_fetch(__sobol_directions, 32*dimension + j);
80         
81         return result;
82 }
83
84 /* lookup index and x/y coordinate, assumes m is a power of two */
85 __device uint sobol_lookup(const uint m, const uint frame, const uint ex, const uint ey, uint *x, uint *y)
86 {
87         /* shift is constant per frame */
88         const uint shift = frame << (m << 1);
89         const uint sobol_shift = sobol(shift);
90         /* van der Corput is its own inverse */
91         const uint lower = van_der_corput(ex << (32 - m));
92         /* need to compensate for ey difference and shift */
93         const uint sobol_lower = sobol(lower);
94         const uint mask = ~-(1 << m) << (32 - m); /* only m upper bits */
95         const uint delta = ((ey << (32 - m)) ^ sobol_lower ^ sobol_shift) & mask;
96         /* only use m upper bits for the index (m is a power of two) */
97         const uint sobol_result = delta | (delta >> m);
98         const uint upper = sobol_inverse(sobol_result);
99         const uint index = shift | upper | lower;
100         *x = van_der_corput(index);
101         *y = sobol_shift ^ sobol_result ^ sobol_lower;
102         return index;
103 }
104
105 __device_inline float path_rng(KernelGlobals *kg, RNG *rng, int sample, int dimension)
106 {
107 #ifdef __SOBOL_FULL_SCREEN__
108         uint result = sobol_dimension(kg, *rng, dimension);
109         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
110         return r;
111 #else
112         /* compute sobol sequence value using direction vectors */
113         uint result = sobol_dimension(kg, sample + SOBOL_SKIP, dimension);
114         float r = (float)result * (1.0f/(float)0xFFFFFFFF);
115
116         /* Cranly-Patterson rotation using rng seed */
117         float shift;
118
119         if(dimension & 1)
120                 shift = (*rng >> 16)*(1.0f/(float)0xFFFF);
121         else
122                 shift = (*rng & 0xFFFF)*(1.0f/(float)0xFFFF);
123
124         return r + shift - floorf(r + shift);
125 #endif
126 }
127
128 __device_inline float path_rng_1D(KernelGlobals *kg, RNG *rng, int sample, int num_samples, int dimension)
129 {
130 #ifdef __CMJ__
131         if(kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_CMJ) {
132                 /* correlated multi-jittered */
133                 int p = *rng + dimension;
134                 return cmj_sample_1D(sample, num_samples, p);
135         }
136 #endif
137
138         /* sobol */
139         return path_rng(kg, rng, sample, dimension);
140 }
141
142 __device_inline float2 path_rng_2D(KernelGlobals *kg, RNG *rng, int sample, int num_samples, int dimension)
143 {
144 #ifdef __CMJ__
145         if(kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_CMJ) {
146                 /* correlated multi-jittered */
147                 int p = *rng + dimension;
148                 return cmj_sample_2D(sample, num_samples, p);
149         }
150 #endif
151
152         /* sobol */
153         return make_float2(path_rng(kg, rng, sample, dimension),
154                            path_rng(kg, rng, sample, dimension + 1));
155 }
156
157 __device_inline void path_rng_init(KernelGlobals *kg, __global uint *rng_state, int sample, int num_samples, RNG *rng, int x, int y, float *fx, float *fy)
158 {
159 #ifdef __SOBOL_FULL_SCREEN__
160         uint px, py;
161         uint bits = 16; /* limits us to 65536x65536 and 65536 samples */
162         uint size = 1 << bits;
163         uint frame = sample;
164
165         *rng = sobol_lookup(bits, frame, x, y, &px, &py);
166
167         *rng ^= kernel_data.integrator.seed;
168
169         if(sample == 0) {
170                 *fx = 0.5f;
171                 *fy = 0.5f;
172         }
173         else {
174                 *fx = size * (float)px * (1.0f/(float)0xFFFFFFFF) - x;
175                 *fy = size * (float)py * (1.0f/(float)0xFFFFFFFF) - y;
176         }
177 #else
178         *rng = *rng_state;
179
180         *rng ^= kernel_data.integrator.seed;
181
182         if(sample == 0) {
183                 *fx = 0.5f;
184                 *fy = 0.5f;
185         }
186         else {
187                 float2 fxy = path_rng_2D(kg, rng, sample, num_samples, PRNG_FILTER_U);
188
189                 *fx = fxy.x;
190                 *fy = fxy.y;
191         }
192 #endif
193 }
194
195 __device void path_rng_end(KernelGlobals *kg, __global uint *rng_state, RNG rng)
196 {
197         /* nothing to do */
198 }
199
200 #else
201
202 /* Linear Congruential Generator */
203
204 __device float path_rng(KernelGlobals *kg, RNG& rng, int sample, int dimension)
205 {
206         /* implicit mod 2^32 */
207         rng = (1103515245*(rng) + 12345);
208         return (float)rng * (1.0f/(float)0xFFFFFFFF);
209 }
210
211 __device_inline float path_rng_1D(KernelGlobals *kg, RNG& rng, int sample, int num_samples, int dimension)
212 {
213         return path_rng(kg, rng, sample, dimension);
214 }
215
216 __device_inline float2 path_rng_2D(KernelGlobals *kg, RNG& rng, int sample, int num_samples, int dimension)
217 {
218         return make_float2(path_rng(kg, rng, sample, dimension),
219                            path_rng(kg, rng, sample, dimension + 1));
220 }
221
222 __device void path_rng_init(KernelGlobals *kg, __global uint *rng_state, int sample, int num_samples, RNG *rng, int x, int y, float *fx, float *fy)
223 {
224         /* load state */
225         *rng = *rng_state;
226
227         *rng ^= kernel_data.integrator.seed;
228
229         if(sample == 0) {
230                 *fx = 0.5f;
231                 *fy = 0.5f;
232         }
233         else {
234                 float2 fxy = path_rng_2D(kg, rng, sample, num_samples, PRNG_FILTER_U);
235
236                 *fx = fxy.x;
237                 *fy = fxy.y;
238         }
239 }
240
241 __device void path_rng_end(KernelGlobals *kg, __global uint *rng_state, RNG rng)
242 {
243         /* store state for next sample */
244         *rng_state = rng;
245 }
246
247 #endif
248
249 __device float lcg_step(uint *rng)
250 {
251         /* implicit mod 2^32 */
252         *rng = (1103515245*(*rng) + 12345);
253         return (float)*rng * (1.0f/(float)0xFFFFFFFF);
254 }
255
256 __device uint lcg_init(uint seed)
257 {
258         uint rng = seed;
259         lcg_step(&rng);
260         return rng;
261 }
262
263 CCL_NAMESPACE_END
264