libutil/random: rewrite whole algorithm

libutil/random.c now includes a custom PRNG, which is the PCG family.
It also overhauls random_uniform(), making it way faster. Names were
changed (s/random/rng32/g), and reentrant versions added.
The PRNG is faster than libc's random().
It is way faster than the previous version, which used

libutil/random.c

@@ -1,46 +1,77 @@
 #include <stdint.h>
-#include <stdlib.h>
+#include <stdio.h>
 #include <time.h>
 
+static uint64_t globalstate;
+
 /*
- * Uniformity is achieved by generating new random numbers until the one
- * returned is outside the range [0, 2**32 % upper_bound).  This
- * guarantees the selected random number will be inside
- * [2**32 % upper_bound, 2**32) which maps back to [0, upper_bound)
- * after reduction modulo upper_bound.
- *
- * Copied off OpenBSD (original is arc4random_uniform)
+ * PCG construction
+ * seeding the RNG means merely setting the initial state.
+ * the increment could also be made part of the seed, just make sure it's odd.
  */
 uint32_t
-random_uniform(uint32_t upper_bound)
+rng32_r(uint64_t *state)
 {
-	uint32_t r, min;
+	uint64_t oldstate = *state;
+	uint32_t r, v;
 
-	if (upper_bound < 2)
-		return 0;
+	*state *= UINT64_C(0x9E3793492EEDC3F7);
+	*state += 0x1337;
 
-	/* 2**32 % x == (2**32 - x) % x */
-	min = -upper_bound % upper_bound;
-
-	/*
-	 * This could theoretically loop forever but each retry has
-	 * p > 0.5 (worst case, usually far better) of selecting a
-	 * number inside the range we need, so it should rarely need
-	 * to re-roll.
-	 */
-	for (;;) {
-		r = random(); /* arc4random() is better, but we don't always have it */
-		if (r >= min)
-			break;
-	}
-
-	return r % upper_bound;
+	r = oldstate >> (64 - 5);
+	v = (oldstate ^ (oldstate >> 18)) >> (32 - 5);
+	v = (v >> (-r & 31)) | (v << r);
+	return v;
 }
 
+uint32_t
+rng32(void)
+{
+	return rng32_r(&globalstate);
+}
+
+/*
+ * Based on optimized Lemire's method
+ * https://pcg-random.org/posts/bounded-rands.html
+ */
+uint32_t
+rng32_bounded_r(uint64_t *state, uint32_t range) {
+	uint32_t x = rng32_r(state);
+	uint64_t m = (uint64_t)x * (uint64_t)range;
+	uint32_t l = (uint32_t)m;
+	if (l < range) {
+		uint32_t t = -range;
+		if (t >= range) {
+			t -= range;
+			if (t >= range)
+				t %= range;
+		}
+		while (l < t) {
+			x = rng32_r(state);
+			m = (uint64_t)x * (uint64_t)range;
+			l = (uint32_t)m;
+		}
+	}
+	return m >> 32;
+}
+
+uint64_t
+rng32_bounded(uint32_t range)
+{
+	return rng32_bounded_r(&globalstate, range);
+}
+
+/* Initialize state with somewhat random number */
 void
-random_seed(void)
+rng32_seed_r(uint64_t *state)
 {
 	struct timespec ts;
 	clock_gettime(CLOCK_REALTIME, &ts);
-	srandom(ts.tv_nsec); /* not a good source of randomness, but eh */
+	*state = (intptr_t)&printf ^ ts.tv_sec ^ ((unsigned long)ts.tv_nsec * 0xAC5533CD);
+}
+
+void
+rng32_seed(void)
+{
+	rng32_seed_r(&globalstate);
 }

util.h

@@ -91,3 +91,6 @@ int mkdirp(const char *, mode_t, mode_t);
 #undef memmem
 #define memmem xmemmem
 void *memmem(const void *, size_t, const void *, size_t);
+uint32_t rng32_r(uint64_t*);
+uint32_t rng32(void);
+uint32_t rng32_bounded_r(uint64_t*, uint32_t);