pcg.cpp 2.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859
  1. // *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
  2. // Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)
  3. #include "pcg.h"
  4. uint32_t pcg32_random_r(pcg32_random_t* rng)
  5. {
  6. uint64_t oldstate = rng->state;
  7. // Advance internal state
  8. rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
  9. // Calculate output function (XSH RR), uses old state for max ILP
  10. uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
  11. uint32_t rot = oldstate >> 59u;
  12. return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
  13. }
  14. // Source from http://www.pcg-random.org/downloads/pcg-c-basic-0.9.zip
  15. void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq)
  16. {
  17. rng->state = 0U;
  18. rng->inc = (initseq << 1u) | 1u;
  19. pcg32_random_r(rng);
  20. rng->state += initstate;
  21. pcg32_random_r(rng);
  22. }
  23. // Source from https://github.com/imneme/pcg-c-basic/blob/master/pcg_basic.c
  24. // pcg32_boundedrand_r(rng, bound):
  25. // Generate a uniformly distributed number, r, where 0 <= r < bound
  26. uint32_t pcg32_boundedrand_r(pcg32_random_t *rng, uint32_t bound) {
  27. // To avoid bias, we need to make the range of the RNG a multiple of
  28. // bound, which we do by dropping output less than a threshold.
  29. // A naive scheme to calculate the threshold would be to do
  30. //
  31. // uint32_t threshold = 0x100000000ull % bound;
  32. //
  33. // but 64-bit div/mod is slower than 32-bit div/mod (especially on
  34. // 32-bit platforms). In essence, we do
  35. //
  36. // uint32_t threshold = (0x100000000ull-bound) % bound;
  37. //
  38. // because this version will calculate the same modulus, but the LHS
  39. // value is less than 2^32.
  40. uint32_t threshold = -bound % bound;
  41. // Uniformity guarantees that this loop will terminate. In practice, it
  42. // should usually terminate quickly; on average (assuming all bounds are
  43. // equally likely), 82.25% of the time, we can expect it to require just
  44. // one iteration. In the worst case, someone passes a bound of 2^31 + 1
  45. // (i.e., 2147483649), which invalidates almost 50% of the range. In
  46. // practice, bounds are typically small and only a tiny amount of the range
  47. // is eliminated.
  48. for (;;) {
  49. uint32_t r = pcg32_random_r(rng);
  50. if (r >= threshold)
  51. return r % bound;
  52. }
  53. }