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