Merge branch 'master' into blender2.8
[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_utildefines.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 static void ensure_rng_thread_safe(void)
272 {
273         /* TODO(sergey): Ideally we will get rid of all rng functions which
274          * are using global generator. But for until then we need some way to
275          * catch "bad" calls at runtime.
276          *
277          * NOTE: Lots of areas are not ported, so we keep check disabled for now.
278          */
279         // BLI_assert(BLI_thread_is_main());
280 }
281
282 float BLI_frand(void)
283 {
284         ensure_rng_thread_safe();
285         return BLI_rng_get_float(&theBLI_rng);
286 }
287
288 float BLI_hash_frand(unsigned int seed)
289 {
290         RNG rng;
291
292         BLI_rng_srandom(&rng, seed);
293         return BLI_rng_get_float(&rng);
294 }
295
296 void BLI_array_randomize(void *data, unsigned int elem_size, unsigned int elem_tot, unsigned int seed)
297 {
298         RNG rng;
299
300         BLI_rng_seed(&rng, seed);
301         BLI_rng_shuffle_array(&rng, data, elem_size, elem_tot);
302 }
303
304 /* ********* for threaded random ************** */
305
306 static RNG rng_tab[BLENDER_MAX_THREADS];
307
308 void BLI_thread_srandom(int thread, unsigned int seed)
309 {
310         if (thread >= BLENDER_MAX_THREADS)
311                 thread = 0;
312
313         BLI_rng_seed(&rng_tab[thread], seed + hash[seed & 255]);
314         seed = BLI_rng_get_uint(&rng_tab[thread]);
315         BLI_rng_seed(&rng_tab[thread], seed + hash[seed & 255]);
316         seed = BLI_rng_get_uint(&rng_tab[thread]);
317         BLI_rng_seed(&rng_tab[thread], seed + hash[seed & 255]);
318 }
319
320 int BLI_thread_rand(int thread)
321 {
322         return BLI_rng_get_int(&rng_tab[thread]);
323 }
324
325 float BLI_thread_frand(int thread)
326 {
327         return BLI_rng_get_float(&rng_tab[thread]);
328 }
329
330 struct RNG_THREAD_ARRAY {
331         RNG rng_tab[BLENDER_MAX_THREADS];
332 };
333
334 RNG_THREAD_ARRAY *BLI_rng_threaded_new(void)
335 {
336         unsigned int i;
337         RNG_THREAD_ARRAY *rngarr = MEM_mallocN(sizeof(RNG_THREAD_ARRAY), "random_array");
338
339         for (i = 0; i < BLENDER_MAX_THREADS; i++) {
340                 BLI_rng_srandom(&rngarr->rng_tab[i], (unsigned int)clock());
341         }
342
343         return rngarr;
344 }
345
346 void BLI_rng_threaded_free(struct RNG_THREAD_ARRAY *rngarr)
347 {
348         MEM_freeN(rngarr);
349 }
350
351 int BLI_rng_thread_rand(RNG_THREAD_ARRAY *rngarr, int thread)
352 {
353         return BLI_rng_get_int(&rngarr->rng_tab[thread]);
354 }
355
356 /* ********* Low-discrepancy sequences ************** */
357
358 /* incremental halton sequence generator, from:
359  * "Instant Radiosity", Keller A. */
360 BLI_INLINE double halton_ex(double invprimes, double *offset)
361 {
362         double e = fabs((1.0 - *offset) - 1e-10);
363
364         if (invprimes >= e) {
365                 double lasth;
366                 double h = invprimes;
367
368                 do {
369                         lasth = h;
370                         h *= invprimes;
371                 } while (h >= e);
372
373                 *offset += ((lasth + h) - 1.0);
374         }
375         else {
376                 *offset += invprimes;
377         }
378
379         return *offset;
380 }
381
382 void BLI_halton_1D(unsigned int prime, double offset, int n, double *r)
383 {
384         const double invprime = 1.0 / (double)prime;
385
386         *r = 0.0;
387
388         for (int s = 0; s < n; s++) {
389                 *r = halton_ex(invprime, &offset);
390         }
391 }
392
393 void BLI_halton_2D(unsigned int prime[2], double offset[2], int n, double *r)
394 {
395         const double invprimes[2] = {1.0 / (double)prime[0], 1.0 / (double)prime[1]};
396
397         r[0] = r[1] = 0.0;
398
399         for (int s = 0; s < n; s++) {
400                 for (int i = 0; i < 2; i++) {
401                         r[i] = halton_ex(invprimes[i], &offset[i]);
402                 }
403         }
404 }
405
406 void BLI_halton_3D(unsigned int prime[3], double offset[3], int n, double *r)
407 {
408         const double invprimes[3] = {1.0 / (double)prime[0], 1.0 / (double)prime[1], 1.0 / (double)prime[2]};
409
410         r[0] = r[1] = r[2] = 0.0;
411
412         for (int s = 0; s < n; s++) {
413                 for (int i = 0; i < 3; i++) {
414                         r[i] = halton_ex(invprimes[i], &offset[i]);
415                 }
416         }
417 }
418
419 void BLI_halton_2D_sequence(unsigned int prime[2], double offset[2], int n, double *r)
420 {
421         const double invprimes[2] = {1.0 / (double)prime[0], 1.0 / (double)prime[1]};
422
423         for (int s = 0; s < n; s++) {
424                 for (int i = 0; i < 2; i++) {
425                         r[s * 2 + i] = halton_ex(invprimes[i], &offset[i]);
426                 }
427         }
428 }
429
430
431 /* From "Sampling with Hammersley and Halton Points" TT Wong
432  * Appendix: Source Code 1 */
433 BLI_INLINE double radical_inverse(unsigned int n)
434 {
435         double u = 0;
436
437         /* This reverse the bitwise representation
438          * around the decimal point. */
439         for (double p = 0.5; n; p *= 0.5, n >>= 1) {
440                 if (n & 1) {
441                         u += p;
442                 }
443         }
444
445         return u;
446 }
447
448 void BLI_hammersley_1D(unsigned int n, double *r)
449 {
450         *r = radical_inverse(n);
451 }
452
453 void BLI_hammersley_2D_sequence(unsigned int n, double *r)
454 {
455         for (unsigned int s = 0; s < n; s++) {
456                 r[s * 2 + 0] = (double)(s + 0.5) / (double)n;
457                 r[s * 2 + 1] = radical_inverse(s);
458         }
459 }