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 }