Cleanup: use BLI_compiler_compat.h for BLI_INLINE
[blender.git] / source / blender / blenlib / intern / rand.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
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  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): none yet.
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/blenlib/intern/rand.c
29  *  \ingroup bli
30  */
31
32
33 #include <stdlib.h>
34 #include <string.h>
35 #include <math.h>
36 #include <time.h>
37
38 #include "MEM_guardedalloc.h"
39
40 #include "BLI_threads.h"
41 #include "BLI_rand.h"
42 #include "BLI_math.h"
43
44 /* defines BLI_INLINE */
45 #include "BLI_compiler_compat.h"
46
47 #include "BLI_sys_types.h"
48 #include "BLI_strict_flags.h"
49
50 #define MULTIPLIER  0x5DEECE66Dll
51 #define MASK        0x0000FFFFFFFFFFFFll
52 #define MASK_BYTES   2
53
54 #define ADDEND      0xB
55 #define LOWSEED     0x330E
56
57 extern unsigned char hash[];    // noise.c
58
59 /**
60  * Random Number Generator.
61  */
62 struct RNG {
63         uint64_t X;
64 };
65
66 RNG *BLI_rng_new(unsigned int seed)
67 {
68         RNG *rng = MEM_mallocN(sizeof(*rng), "rng");
69
70         BLI_rng_seed(rng, seed);
71
72         return rng;
73 }
74
75 /**
76  * A version of #BLI_rng_new that hashes the seed.
77  */
78 RNG *BLI_rng_new_srandom(unsigned int seed)
79 {
80         RNG *rng = MEM_mallocN(sizeof(*rng), "rng");
81
82         BLI_rng_srandom(rng, seed);
83
84         return rng;
85 }
86
87 void BLI_rng_free(RNG *rng)
88 {
89         MEM_freeN(rng);
90 }
91
92 void BLI_rng_seed(RNG *rng, unsigned int seed)
93 {
94         rng->X = (((uint64_t) seed) << 16) | LOWSEED;
95 }
96
97 /**
98  * Use a hash table to create better seed.
99  */
100 void BLI_rng_srandom(RNG *rng, unsigned int seed)
101 {
102         BLI_rng_seed(rng, seed + hash[seed & 255]);
103         seed = BLI_rng_get_uint(rng);
104         BLI_rng_seed(rng, seed + hash[seed & 255]);
105         seed = BLI_rng_get_uint(rng);
106         BLI_rng_seed(rng, seed + hash[seed & 255]);
107 }
108
109 BLI_INLINE void rng_step(RNG *rng)
110 {
111         rng->X = (MULTIPLIER * rng->X + ADDEND) & MASK;
112 }
113
114 void BLI_rng_get_char_n(RNG *rng, char *bytes, size_t bytes_len)
115 {
116         size_t last_len = 0;
117         size_t trim_len = bytes_len;
118
119 #define RAND_STRIDE (sizeof(rng->X) - MASK_BYTES)
120
121         if (trim_len > RAND_STRIDE) {
122                 last_len = trim_len % RAND_STRIDE;
123                 trim_len = trim_len - last_len;
124         }
125         else {
126                 trim_len = 0;
127                 last_len = bytes_len;
128         }
129
130         const char *data_src = (void *)&(rng->X);
131         size_t i = 0;
132         while (i != trim_len) {
133                 BLI_assert(i < trim_len);
134 #ifdef __BIG_ENDIAN__
135                 for (size_t j = (RAND_STRIDE + MASK_BYTES) - 1; j != MASK_BYTES - 1; j--)
136 #else
137                 for (size_t j = 0; j != RAND_STRIDE; j++)
138 #endif
139                 {
140                         bytes[i++] = data_src[j];
141                 }
142                 rng_step(rng);
143         }
144         if (last_len) {
145                 for (size_t j = 0; j != last_len; j++) {
146                         bytes[i++] = data_src[j];
147                 }
148         }
149
150 #undef RAND_STRIDE
151 }
152
153 int BLI_rng_get_int(RNG *rng)
154 {
155         rng_step(rng);
156         return (int) (rng->X >> 17);
157 }
158
159 unsigned int BLI_rng_get_uint(RNG *rng)
160 {
161         rng_step(rng);
162         return (unsigned int) (rng->X >> 17);
163 }
164
165 /**
166  * \return Random value (0..1), but never 1.0.
167  */
168 double BLI_rng_get_double(RNG *rng)
169 {
170         return (double) BLI_rng_get_int(rng) / 0x80000000;
171 }
172
173 /**
174  * \return Random value (0..1), but never 1.0.
175  */
176 float BLI_rng_get_float(RNG *rng)
177 {
178         return (float) BLI_rng_get_int(rng) / 0x80000000;
179 }
180
181 void BLI_rng_get_float_unit_v2(RNG *rng, float v[2])
182 {
183         float a = (float)(M_PI * 2.0) * BLI_rng_get_float(rng);
184         v[0] = cosf(a);
185         v[1] = sinf(a);
186 }
187
188 void BLI_rng_get_float_unit_v3(RNG *rng, float v[3])
189 {
190         float r;
191         v[2] = (2.0f * BLI_rng_get_float(rng)) - 1.0f;
192         if ((r = 1.0f - (v[2] * v[2])) > 0.0f) {
193                 float a = (float)(M_PI * 2.0) * BLI_rng_get_float(rng);
194                 r = sqrtf(r);
195                 v[0] = r * cosf(a);
196                 v[1] = r * sinf(a);
197         }
198         else {
199                 v[2] = 1.0f;
200         }
201 }
202
203 /**
204  * Generate a random point inside given tri.
205  */
206 void BLI_rng_get_tri_sample_float_v2(
207         RNG *rng, const float v1[2], const float v2[2], const float v3[2],
208         float r_pt[2])
209 {
210         float u = BLI_rng_get_float(rng);
211         float v = BLI_rng_get_float(rng);
212
213         float side_u[2], side_v[2];
214
215         if ((u + v) > 1.0f) {
216                 u = 1.0f - u;
217                 v = 1.0f - v;
218         }
219
220         sub_v2_v2v2(side_u, v2, v1);
221         sub_v2_v2v2(side_v, v3, v1);
222
223         copy_v2_v2(r_pt, v1);
224         madd_v2_v2fl(r_pt, side_u, u);
225         madd_v2_v2fl(r_pt, side_v, v);
226 }
227
228 void BLI_rng_shuffle_array(RNG *rng, void *data, unsigned int elem_size_i, unsigned int elem_tot)
229 {
230         const size_t elem_size = (size_t)elem_size_i;
231         unsigned int i = elem_tot;
232         void *temp;
233
234         if (elem_tot <= 1) {
235                 return;
236         }
237
238         temp = malloc(elem_size);
239
240         while (i--) {
241                 unsigned int j = BLI_rng_get_uint(rng) % elem_tot;
242                 if (i != j) {
243                         void *iElem = (unsigned char *)data + i * elem_size_i;
244                         void *jElem = (unsigned char *)data + j * elem_size_i;
245                         memcpy(temp, iElem, elem_size);
246                         memcpy(iElem, jElem, elem_size);
247                         memcpy(jElem, temp, elem_size);
248                 }
249         }
250
251         free(temp);
252 }
253
254 /**
255  * Simulate getting \a n random values.
256  *
257  * \note Useful when threaded code needs consistent values, independent of task division.
258  */
259 void BLI_rng_skip(RNG *rng, int n)
260 {
261         while (n--) {
262                 rng_step(rng);
263         }
264 }
265
266 /***/
267
268 /* initialize with some non-zero seed */
269 static RNG theBLI_rng = {611330372042337130};
270
271 void BLI_srandom(unsigned int seed)
272 {
273         BLI_rng_srandom(&theBLI_rng, seed);
274 }
275
276 int BLI_rand(void)
277 {
278         return BLI_rng_get_int(&theBLI_rng);
279 }
280
281 float BLI_frand(void)
282 {
283         return BLI_rng_get_float(&theBLI_rng);
284 }
285
286 void BLI_frand_unit_v3(float v[3])
287 {
288         BLI_rng_get_float_unit_v3(&theBLI_rng, v);
289 }
290
291 float BLI_hash_frand(unsigned int seed)
292 {
293         RNG rng;
294
295         BLI_rng_srandom(&rng, seed);
296         return BLI_rng_get_float(&rng);
297 }
298
299 void BLI_array_randomize(void *data, unsigned int elem_size, unsigned int elem_tot, unsigned int seed)
300 {
301         RNG rng;
302
303         BLI_rng_seed(&rng, seed);
304         BLI_rng_shuffle_array(&rng, data, elem_size, elem_tot);
305 }
306
307 /* ********* for threaded random ************** */
308
309 static RNG rng_tab[BLENDER_MAX_THREADS];
310
311 void BLI_thread_srandom(int thread, unsigned int seed)
312 {
313         if (thread >= BLENDER_MAX_THREADS)
314                 thread = 0;
315
316         BLI_rng_seed(&rng_tab[thread], seed + hash[seed & 255]);
317         seed = BLI_rng_get_uint(&rng_tab[thread]);
318         BLI_rng_seed(&rng_tab[thread], seed + hash[seed & 255]);
319         seed = BLI_rng_get_uint(&rng_tab[thread]);
320         BLI_rng_seed(&rng_tab[thread], seed + hash[seed & 255]);
321 }
322
323 int BLI_thread_rand(int thread)
324 {
325         return BLI_rng_get_int(&rng_tab[thread]);
326 }
327
328 float BLI_thread_frand(int thread)
329 {
330         return BLI_rng_get_float(&rng_tab[thread]);
331 }
332
333 struct RNG_THREAD_ARRAY {
334         RNG rng_tab[BLENDER_MAX_THREADS];
335 };
336
337 RNG_THREAD_ARRAY *BLI_rng_threaded_new(void)
338 {
339         unsigned int i;
340         RNG_THREAD_ARRAY *rngarr = MEM_mallocN(sizeof(RNG_THREAD_ARRAY), "random_array");
341
342         for (i = 0; i < BLENDER_MAX_THREADS; i++) {
343                 BLI_rng_srandom(&rngarr->rng_tab[i], (unsigned int)clock());
344         }
345
346         return rngarr;
347 }
348
349 void BLI_rng_threaded_free(struct RNG_THREAD_ARRAY *rngarr)
350 {
351         MEM_freeN(rngarr);
352 }
353
354 int BLI_rng_thread_rand(RNG_THREAD_ARRAY *rngarr, int thread)
355 {
356         return BLI_rng_get_int(&rngarr->rng_tab[thread]);
357 }
358
359 /* ********* Low-discrepancy sequences ************** */
360
361 /* incremental halton sequence generator, from:
362  * "Instant Radiosity", Keller A. */
363 BLI_INLINE double halton_ex(double invprimes, double *offset)
364 {
365         double e = fabs((1.0 - *offset) - 1e-10);
366
367         if (invprimes >= e) {
368                 double lasth;
369                 double h = invprimes;
370
371                 do {
372                         lasth = h;
373                         h *= invprimes;
374                 } while (h >= e);
375
376                 *offset += ((lasth + h) - 1.0);
377         }
378         else {
379                 *offset += invprimes;
380         }
381
382         return *offset;
383 }
384
385 void BLI_halton_1D(unsigned int prime, double offset, int n, double *r)
386 {
387         const double invprime = 1.0 / (double)prime;
388
389         for (int s = 0; s < n; s++) {
390                 *r = halton_ex(invprime, &offset);
391         }
392 }
393
394 void BLI_halton_2D(unsigned int prime[2], double offset[2], int n, double *r)
395 {
396         const double invprimes[2] = {1.0 / (double)prime[0], 1.0 / (double)prime[1]};
397
398         for (int s = 0; s < n; s++) {
399                 for (int i = 0; i < 2; i++) {
400                         r[i] = halton_ex(invprimes[i], &offset[i]);
401                 }
402         }
403 }
404
405 void BLI_halton_3D(unsigned int prime[3], double offset[3], int n, double *r)
406 {
407         const double invprimes[3] = {1.0 / (double)prime[0], 1.0 / (double)prime[1], 1.0 / (double)prime[2]};
408
409         for (int s = 0; s < n; s++) {
410                 for (int i = 0; i < 3; i++) {
411                         r[i] = halton_ex(invprimes[i], &offset[i]);
412                 }
413         }
414 }
415
416 void BLI_halton_2D_sequence(unsigned int prime[2], double offset[2], int n, double *r)
417 {
418         const double invprimes[2] = {1.0 / (double)prime[0], 1.0 / (double)prime[1]};
419
420         for (int s = 0; s < n; s++) {
421                 for (int i = 0; i < 2; i++) {
422                         r[s * 2 + i] = halton_ex(invprimes[i], &offset[i]);
423                 }
424         }
425 }
426
427
428 /* From "Sampling with Hammersley and Halton Points" TT Wong
429  * Appendix: Source Code 1 */
430 BLI_INLINE double radical_inverse(unsigned int n)
431 {
432         double u = 0;
433
434         /* This reverse the bitwise representation
435          * around the decimal point. */
436         for (double p = 0.5; n; p *= 0.5, n >>= 1) {
437                 if (n & 1) {
438                         u += p;
439                 }
440         }
441
442         return u;
443 }
444
445 void BLI_hammersley_1D(unsigned int n, double *r)
446 {
447         *r = radical_inverse(n);
448 }
449
450 void BLI_hammersley_2D_sequence(unsigned int n, double *r)
451 {
452         for (unsigned int s = 0; s < n; s++) {
453                 r[s * 2 + 0] = (double)(s + 0.5) / (double)n;
454                 r[s * 2 + 1] = radical_inverse(s);
455         }
456 }