2 * UCW Library -- Fast pseudo-random generator
4 * (c) 2020--2022 Pavel Charvat <pchar@ucw.cz>
6 * This software may be freely distributed and used according to the terms
7 * of the GNU Lesser General Public License.
13 #include <ucw/random.h>
14 #include <ucw/threads.h>
15 #include <ucw/unaligned.h>
24 // Select exactly one of these pseudo-random generators:
25 //#define FASTRAND_USE_PCG
26 //#define FASTRAND_USE_SPLITMIX
27 #define FASTRAND_USE_XOROSHIRO
28 //#define FASTRAND_USE_LIBC_SVID
34 * -- https://en.wikipedia.org/wiki/Permuted_congruential_generator
35 * -- https://www.pcg-random.org/
36 * -- "PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation",
37 * 2014, Melissa E. O’Neill
40 * -- Simple and fast generator with good statistical properties.
41 * -- We use a constant increment and 64->32bit XSH RR output function.
44 #define FASTRAND_PCG_MUL 6364136223846793005LLU
45 #define FASTRAND_PCG_INC 1442695040888963407LLU
48 * FASTRAND_USE_SPLITMIX
51 * -- http://prng.di.unimi.it/splitmix64.c (public, no copyright)
52 * -- http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
53 * -- "Fast Splittable Pseudorandom Number Generators",
54 * 2014, Guy L. Steele Jr., Doug Lea, Christine H. Flood
57 * -- Each 64 bit number is generated exactly once in a cycle,
58 * so we only use half of it -> 32 bit output.
62 * FASTRAND_USE_XOROSHIRO
65 * -- http://prng.di.unimi.it/ (examples are public, no copyright)
66 * -- "Scrambled Linear Pseudorandom Number Generators",
67 * 2019, David Blackman, Sebastiano Vigna
70 * -- A family of very fast generators. Good statistical properties,
71 * some variants just have problems with lowest bits.
72 * -- We support more xo(ro)shiro variants, see the options below.
75 // Engine: xoroshiro 128, xoshiro 256
76 #define FASTRAND_XOROSHIRO_SIZE 128
78 // Scrambler: 1 (+), 2 (++), 3 (**)
79 #define FASTRAND_XOROSHIRO_SCRAMBLER 1
81 // If defined, we only use the highest 60 bits from each 64 bit number;
82 // useful especially for the weakest '+' scrambler. Even if not
83 // we prefer higher bits where trivial.
84 #undef FASTRAND_XOROSHIRO_ONLY_HIGHER
87 * FASTRAND_USE_LIBC_SVID
89 * -- libc SVID generator
93 #ifdef FASTRAND_USE_PCG
94 u64 pcg_state; // PCG generator's current state
97 #ifdef FASTRAND_USE_SPLITMIX
101 #ifdef FASTRAND_USE_XOROSHIRO
102 u64 xoroshiro_state[FASTRAND_XOROSHIRO_SIZE / 64];
105 #ifdef FASTRAND_USE_LIBC_SVID
106 struct drand48_data svid_data; // SVID generator's state structure
109 bool inited; // Is seeded?
110 fastrand_seed_t seed_value; // Initial seed value; useful for reproducible debugging
113 /*** Common helpers ***/
115 // GCC should optimize these to a single instruction
117 static inline u32 fastrand_ror32(u32 x, uint r)
119 return (x >> r) | (x << (-r & 31));
122 static inline u64 fastrand_rol64(u64 x, uint r)
124 return (x << r) | (x >> (-r & 63));
127 /*** PCG random generator ***/
129 #ifdef FASTRAND_USE_PCG
130 #define FASTRAND_ONE_BITS 32
132 static inline u32 fastrand_one(struct fastrand *c)
134 u64 s = c->pcg_state;
136 c->pcg_state = s * FASTRAND_PCG_MUL + FASTRAND_PCG_INC;
137 u32 x = ((s >> 18) ^ s) >> 27;
138 return fastrand_ror32(x, r);
141 static void fastrand_init(struct fastrand *c)
143 c->pcg_state = c->seed_value + FASTRAND_PCG_INC;
144 (void) fastrand_one(c);
145 DBG("RANDOM: Initialized PCG random generator, seed=0x%jx", (uintmax_t)c->seed_value);
149 /*** splitmix64 generator ***/
151 // A helper; used also for seeding xoroshiro
152 static inline u64 fastrand_splitmix64_next(u64 *state)
154 u64 z = (*state += 0x9e3779b97f4a7c15LLU);
155 z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9LLU;
156 z = (z ^ (z >> 27)) * 0x94d049bb133111ebLLU;
157 return z ^ (z >> 31);
160 #ifdef FASTRAND_USE_SPLITMIX
161 #define FASTRAND_ONE_BITS 32
163 static void fastrand_init(struct fastrand *c)
165 c->splitmix_state = c->seed_value;
166 DBG("RANDOM: Initialized splitmix64 random generator, seed=0x%jx", (uintmax_t)c->seed_value);
169 static inline u32 fastrand_one(struct fastrand *c)
171 return fastrand_splitmix64_next(&c->splitmix_state);
175 /*** xoroshiro variants ***/
177 #ifdef FASTRAND_USE_XOROSHIRO
178 #ifdef FASTRAND_XOROSHIRO_ONLY_HIGHER
179 # define FASTRAND_ONE_BITS 60
181 # define FASTRAND_ONE_BITS 64
182 # define FASTRAND_PREFER_HIGHER
185 static void fastrand_init(struct fastrand *c)
187 u64 x = c->seed_value;
188 for (uint i = 0; i < ARRAY_SIZE(c->xoroshiro_state); i++)
189 c->xoroshiro_state[i] = fastrand_splitmix64_next(&x);
190 if (unlikely(!c->xoroshiro_state[0]))
191 c->xoroshiro_state[0] = 1;
192 DBG("RANDOM: Initialized xoroshiro random generator, seed=0x%jx, size=%u, scrambler=%u",
193 (uintmax_t)c->seed_value, FASTRAND_XOROSHIRO_SIZE, FASTRAND_XOROSHIRO_SCRAMBLER);
196 static inline u64 fastrand_one(struct fastrand *c)
199 u64 *s = c->xoroshiro_state;
201 #if FASTRAND_XOROSHIRO_SIZE == 128
204 #if FASTRAND_XOROSHIRO_SCRAMBLER == 1
206 uint ra = 24, rb = 16, rc = 37;
207 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 2
208 result = fastrand_rol64(s0 + s1, 17) + s0;
209 uint ra = 49, rb = 21, rc = 28;
210 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 3
211 result = fastrand_rol64(s0 * 5, 7) * 9;
212 uint ra = 24, rb = 16, rc = 37;
214 # error "Unsupported xoroshiro parameters"
216 s[0] = fastrand_rol64(s0, ra) ^ s1 ^ (s1 << rb);
217 s[1] = fastrand_rol64(s1, rc);
219 #elif FASTRAND_XOROSHIRO_SIZE == 256
220 #if FASTRAND_XOROSHIRO_SCRAMBLER == 1
221 result = s[0] + s[3];
222 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 2
223 result = fastrand_rol64(s[0] + s[3], 23) + s[0];
224 #elif FASTRAND_XOROSHIRO_SCRAMBLER == 3
225 result = fastrand_rol64(s[1] * 5, 7) * 9;
227 # error "Unsupported xoroshiro parameters"
235 s[3] = fastrand_rol64(s[3], 45);
237 # error "Unsupported xoroshiro parameters"
240 #ifdef FASTRAND_XOROSHIRO_ONLY_HIGHER
248 /*** Libc SVID generator ***/
250 #ifdef FASTRAND_USE_LIBC_SVID
251 #define FASTRAND_ONE_BITS 32
253 static void fastrand_init(struct fastrand *c)
255 DBG("RANDOM: Initialized libc SVID random generator, seed=0x%jx", (uintmax_t)c->seed_value);
256 if (srand48_r(c->seed_value, &c->svid_data) < 0)
257 die("RANDOM: Failed to call srand48_r(): %m");
261 static NONRET void fastrand_mrand48_r_failed(void)
263 die("RANDOM: Failed to call mrand48_r(): %m");
267 static inline u32 fastrand_one(struct fastrand *c)
271 if (unlikely(mrand48_r(&c->svid_data, &r) < 0))
272 fastrand_mrand48_r_failed();
274 (void) mrand48_r(&c->svid_data, &r);
281 /*** Interface for creating/deleting/setting contexts ***/
283 struct fastrand *fastrand_new(void)
285 struct fastrand *c = xmalloc_zero(sizeof(*c));
289 void fastrand_delete(struct fastrand *c)
295 void fastrand_set_seed(struct fastrand *c, fastrand_seed_t seed)
297 c->seed_value = seed;
302 fastrand_seed_t fastrand_gen_seed_value(void)
304 // Generate a seed value as a mixture of:
306 // -- process id to deal with frequent fork() or startups
307 // -- counter to deal with multiple contexts, for example in threaded application
309 static u64 static_counter;
311 val = (static_counter += 0x1927a7639274162dLLU);
315 if (gettimeofday(&tv, NULL) < 0)
316 die("RANDOM: Failed to call gettimeofday(): %m");
317 val += (u64)(tv.tv_sec * 1000000LL + tv.tv_usec) * 0x5a03173d83f78347LLU;
321 fastrand_seed_t fastrand_gen_seed(struct fastrand *c)
323 fastrand_seed_t seed = fastrand_gen_seed_value();
324 fastrand_set_seed(c, seed);
328 bool fastrand_has_seed(struct fastrand *c)
333 void fastrand_reset(struct fastrand *c)
339 // XXX: Adds some overhead, but probably worth it
340 #define FASTRAND_CHECK_SEED(c) do { if (unlikely(!(c)->inited)) fastrand_check_seed_failed(); } while (0)
341 static NONRET void fastrand_check_seed_failed(void)
346 #define FASTRAND_CHECK_SEED(c) do { } while (0)
349 /*** Reading interface ***/
351 // We assume 32-64 bit output of fastrand_one().
352 // For 32 bits fastrand_one() can return u32 type, otherwise it must be u64.
353 COMPILE_ASSERT(fastrand_one_range, FASTRAND_ONE_BITS >= 32 && FASTRAND_ONE_BITS <= 64);
355 #if FASTRAND_ONE_BITS == 32
356 #define FASTRAND_ONE_MAX UINT32_MAX
358 #define FASTRAND_ONE_MAX (UINT64_MAX >> (64 - FASTRAND_ONE_BITS))
361 static inline u32 fastrand_u32_helper(struct fastrand *c)
363 #if FASTRAND_ONE_BITS > 32 && defined(FASTRAND_PREFER_HIGHER)
364 return fastrand_one(c) >> (FASTRAND_ONE_BITS - 32);
366 return fastrand_one(c);
370 u32 fastrand_u32(struct fastrand *c)
372 FASTRAND_CHECK_SEED(c);
373 return fastrand_u32_helper(c);
376 u64 fastrand_u64(struct fastrand *c)
378 FASTRAND_CHECK_SEED(c);
379 #if FASTRAND_ONE_BITS < 64
380 return ((u64)fastrand_one(c) << FASTRAND_ONE_BITS) | fastrand_one(c);
382 return fastrand_one(c);
386 u64 fastrand_bits_u64(struct fastrand *c, uint bits)
391 // Notice that we never rotate an integer by all its bits (forbidden by C standard).
392 FASTRAND_CHECK_SEED(c);
393 #if FASTRAND_ONE_BITS == 64
396 return fastrand_one(c) >> (FASTRAND_ONE_BITS - bits);
398 u64 r = fastrand_one(c);
399 if (bits <= FASTRAND_ONE_BITS)
400 return r >> (FASTRAND_ONE_BITS - bits);
401 return (r << (bits - FASTRAND_ONE_BITS)) | (fastrand_one(c) >> (2 * FASTRAND_ONE_BITS - bits));
406 * fastrand_max() and variants
410 * We can return r % bound from uniformly generated i-bit number r,
412 * r - r % bound + bound <= 2^i
413 * r - r % bound <= 2^i-1 - bound + 1 -- no risk of integer overflow
414 * r - r % bound <= (2^i-1 + 1) - bound -- can temporarily overflow
416 * Then r fits to a complete interval in the i-bit range
417 * and results should be uniformly distributed.
419 * If r is too high, we simply generate a different r and repeat
420 * the whole process. The probability is < 50% for any bound <= 2^i
421 * (with worst case bound == 2^i / 2 + 1), so the average number
422 * of iterations is < 2.
425 * Algorithm 2 (without any division in many cases, but needs big multiplication):
426 * ("Fast Random Integer Generation in an Interval", 2019, Daniel Lemier)
428 * Let bound < 2^i. Then for j in [0, bound) each interval [j * 2^i + 2^i % bound, (j + 1) * 2^i)
429 * contains the same number of multiples [0, 1 * bound, 2 * bound, ..., 2^i * bound),
430 * exactly 2^i / bound (rounded down) of them.
432 * Therefore if we have a random i-bit number r and r * bound fits to one of these intervals,
433 * we can immediately return (r * bound) / 2^i (index of that interval) and results should be
434 * correctly distributed. If we are lucky and r fits to even smaller interval
435 * [j * 2^i + bound, (j + 1) * 2^i], we can check that very quickly without any division.
437 * Probability of multiple iterations is again < 50%.
438 * Probability of division-less check is (2^i - bound) / 2^i.
441 * Algorithm 3 (always without divisions or multiplications, but needs more iterations):
443 * Let 2^(i-1) <= bound < 2^i. Then i-bit random number r has at least 50% probability
444 * to be lower than bound, when we can simply return it.
445 * Can be useful for very fast fastrand_one() like xoroshiro128+.
448 #define FASTRAND_MAX_U32_ALGORITHM 2
449 #define FASTRAND_MAX_U64_ALGORITHM 3
451 static inline u32 fastrand_max_u32_helper(struct fastrand *c, u32 bound)
453 // Variant for bound <= UINT32_MAX.
454 // This is the most common case -> we optimize it and inline the code to both fastrand_max() and fastrand_max_u64().
456 #if FASTRAND_MAX_U32_ALGORITHM == 1
459 u32 r = fastrand_u32_helper(c);
461 if (r - x <= (u32)-bound)
464 // Possible optimization for unlucky cases
467 r = fastrand_u32_helper(c);
473 #elif FASTRAND_MAX_U32_ALGORITHM == 2
474 u32 r = fastrand_u32_helper(c);
475 u64 m = (u64)r * bound;
481 // This condition can decrease the probability of division to <= 50% even in worst case and also optimize bad cases.
484 // Implies l < t % bound, so we don't need to check it before switching to completely different loop.
486 r = fastrand_u32_helper(c);
494 r = fastrand_u32_helper(c);
501 #elif FASTRAND_MAX_U32_ALGORITHM == 3
502 // XXX: Probably could be optimized with __builtin_clz()
503 u32 r, y = bound - 1;
510 r = fastrand_u32_helper(c) & y;
515 # error Invalid FASTRAND_MAX_U32_ALGORITHM
519 static u64 fastrand_max_u64_helper(struct fastrand *c, u64 bound)
521 // Less common case -> compromise between speed and size
522 #if FASTRAND_MAX_U64_ALGORITHM == 1
525 u64 r = fastrand_one(c);
526 u64 m = FASTRAND_ONE_MAX;
527 #if FASTRAND_ONE_BITS < 64
530 r = (r << FASTRAND_ONE_BITS) | fastrand_one(c);
531 m = (m << FASTRAND_ONE_BITS) | FASTRAND_ONE_MAX;
535 if (r - x <= m + 1 - bound)
539 #elif FASTRAND_MAX_U64_ALGORITHM == 3
540 u64 r, y = bound - 1;
550 #if FASTRAND_ONE_BITS < 64
551 if (y > FASTRAND_ONE_MAX)
552 r = (r << FASTRAND_ONE_BITS) | fastrand_one(c);
559 # error Invalid FASTRAND_MAX_U64_ALGORITHM
563 uint fastrand_max(struct fastrand *c, uint bound)
565 FASTRAND_CHECK_SEED(c);
566 #if UINT32_MAX < UINT_MAX
567 if (bound <= UINT32_MAX)
568 return fastrand_max_u32_helper(c, bound);
570 return fastrand_max_u64_helper(c, bound);
572 return fastrand_max_u32_helper(c, bound);
576 u64 fastrand_max_u64(struct fastrand *c, u64 bound)
578 FASTRAND_CHECK_SEED(c);
579 if (bound <= UINT32_MAX)
580 return fastrand_max_u32_helper(c, bound);
582 return fastrand_max_u64_helper(c, bound);
585 void fastrand_mem(struct fastrand *c, void *ptr, size_t size)
587 // XXX: Could be optimized to aligned assignments instead of put_u32/64, but it
588 // would be quite complex to make it deterministic (for same seed and any @ptr)
589 // and large sizes are not used often anyway.
590 FASTRAND_CHECK_SEED(c);
592 #if FASTRAND_ONE_BITS < 64
593 for (size_t n = size >> 2; n--; p += 4)
594 put_u32(p, fastrand_u32_helper(c));
597 u32 x = fastrand_u32_helper(c);
599 for (size_t n = size >> 3; n--; p += 8)
600 put_u64(p, fastrand_one(c));
603 u64 x = fastrand_one(c);
613 double fastrand_double(struct fastrand *c)
615 // Generate 52 random bits
616 FASTRAND_CHECK_SEED(c);
617 u64 r = fastrand_one(c);
618 #if FASTRAND_ONE_BITS >= 52
619 r >>= FASTRAND_ONE_BITS - 52;
621 // Do we want to generate doubles with full 52 precision?
623 r = (r << (52 - FASTRAND_ONE_BITS)) | (fastrand_one(c) >> (2 * FASTRAND_ONE_BITS - 52));
625 r <<= 52 - FASTRAND_ONE_BITS;
629 // With standard IEC 60559 / IEEE 754 format, we can directly store @r as a double
630 // from range [1.0, 2.0) and subtract 1.0.
631 // But the alternative with multiplication is often faster anyway.
632 #if 0 && defined(__STDC_IEC_559__)
637 ASSERT(sizeof(double) == 8 && sizeof(tmp) == 8);
638 tmp.u = (0x3ffLLU << 52) | r;
641 static const double mul = 1.0 / (1LLU << 52);
644 return ldexp(r, -52);