pcg_basic.c 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117
  1. /*
  2. * PCG Random Number Generation for C.
  3. *
  4. * Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
  5. *
  6. * Licensed under the Apache License, Version 2.0 (the "License");
  7. * you may not use this file except in compliance with the License.
  8. * You may obtain a copy of the License at
  9. *
  10. * http://www.apache.org/licenses/LICENSE-2.0
  11. *
  12. * Unless required by applicable law or agreed to in writing, software
  13. * distributed under the License is distributed on an "AS IS" BASIS,
  14. * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  15. * See the License for the specific language governing permissions and
  16. * limitations under the License.
  17. *
  18. * For additional information about the PCG random number generation scheme,
  19. * including its license and other licensing options, visit
  20. *
  21. * http://www.pcg-random.org
  22. */
  23. /*
  24. * This code is derived from the full C implementation, which is in turn
  25. * derived from the canonical C++ PCG implementation. The C++ version
  26. * has many additional features and is preferable if you can use C++ in
  27. * your project.
  28. */
  29. #include "pcg_basic.h"
  30. // state for global RNGs
  31. static pcg32_random_t pcg32_global = PCG32_INITIALIZER;
  32. // pcg32_srandom(initial_state, initseq)
  33. // pcg32_srandom_r(rng, initial_state, initseq):
  34. // Seed the rng. Specified in two parts, state initializer and a
  35. // sequence selection constant (a.k.a. stream id)
  36. void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initial_state, uint64_t initseq)
  37. {
  38. rng->state = 0U;
  39. rng->inc = (initseq << 1u) | 1u;
  40. pcg32_random_r(rng);
  41. rng->state += initial_state;
  42. pcg32_random_r(rng);
  43. }
  44. void pcg32_srandom(uint64_t seed, uint64_t seq)
  45. {
  46. pcg32_srandom_r(&pcg32_global, seed, seq);
  47. }
  48. // pcg32_random()
  49. // pcg32_random_r(rng)
  50. // Generate a uniformly distributed 32-bit random number
  51. uint32_t pcg32_random_r(pcg32_random_t* rng)
  52. {
  53. uint64_t oldstate = rng->state;
  54. rng->state = oldstate * 6364136223846793005ULL + rng->inc;
  55. uint32_t xorshifted = (uint32_t)(((oldstate >> 18u) ^ oldstate) >> 27u);
  56. uint32_t rot = (uint32_t)(oldstate >> 59u);
  57. return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
  58. }
  59. uint32_t pcg32_random()
  60. {
  61. return pcg32_random_r(&pcg32_global);
  62. }
  63. // pcg32_boundedrand(bound):
  64. // pcg32_boundedrand_r(rng, bound):
  65. // Generate a uniformly distributed number, r, where 0 <= r < bound
  66. uint32_t pcg32_boundedrand_r(pcg32_random_t* rng, uint32_t bound)
  67. {
  68. // To avoid bias, we need to make the range of the RNG a multiple of
  69. // bound, which we do by dropping output less than a threshold.
  70. // A naive scheme to calculate the threshold would be to do
  71. //
  72. // uint32_t threshold = 0x100000000ull % bound;
  73. //
  74. // but 64-bit div/mod is slower than 32-bit div/mod (especially on
  75. // 32-bit platforms). In essence, we do
  76. //
  77. // uint32_t threshold = (0x100000000ull-bound) % bound;
  78. //
  79. // because this version will calculate the same modulus, but the LHS
  80. // value is less than 2^32.
  81. uint32_t threshold = -bound % bound;
  82. // Uniformity guarantees that this loop will terminate. In practice, it
  83. // should usually terminate quickly; on average (assuming all bounds are
  84. // equally likely), 82.25% of the time, we can expect it to require just
  85. // one iteration. In the worst case, someone passes a bound of 2^31 + 1
  86. // (i.e., 2147483649), which invalidates almost 50% of the range. In
  87. // practice, bounds are typically small and only a tiny amount of the range
  88. // is eliminated.
  89. for (;;) {
  90. uint32_t r = pcg32_random_r(rng);
  91. if (r >= threshold)
  92. return r % bound;
  93. }
  94. }
  95. uint32_t pcg32_boundedrand(uint32_t bound)
  96. {
  97. return pcg32_boundedrand_r(&pcg32_global, bound);
  98. }