Replace C++ RNG with xoroshiro128+ and Box-Muller

include/util.hpp

@@ -58,7 +58,9 @@ T *construct(F f, V v, Args... args) {
 // RNG
 ////////////////////
 
+void randomSeedTime();
 uint32_t randomu32();
+uint64_t randomu64();
 /** Returns a uniform random float in the interval [0.0, 1.0) */
 float randomf();
 /** Returns a normal random number with mean 0 and std dev 1 */

src/main.cpp

@@ -10,6 +10,8 @@
 using namespace rack;
 
 int main(int argc, char* argv[]) {
+	randomSeedTime();
+
 	{
 		char *cwd = getcwd(NULL, 0);
 		printf("Current working directory: %s\n", cwd);

src/util.cpp

@@ -4,6 +4,7 @@
 #include <string.h>
 #include <random>
 #include <libgen.h> // for dirname and basename
+#include <sys/time.h>
 
 #if ARCH_WIN
 #include <windows.h>
@@ -12,24 +13,70 @@
 
 namespace rack {
 
-// TODO
-// Convert this to xoroshiro128+ or something, and write custom normal dist implementation
+
+// xoroshiro128+
+// from http://xoroshiro.di.unimi.it/xoroshiro128plus.c
+
+static uint64_t xoroshiro128plus_state[2];
+
+static uint64_t rotl(const uint64_t x, int k) {
+	return (x << k) | (x >> (64 - k));
+}
+
+static uint64_t xoroshiro128plus_next(void) {
+	const uint64_t s0 = xoroshiro128plus_state[0];
+	uint64_t s1 = xoroshiro128plus_state[1];
+	const uint64_t result = s0 + s1;
+
+	s1 ^= s0;
+	xoroshiro128plus_state[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); // a, b
+	xoroshiro128plus_state[1] = rotl(s1, 36); // c
+
+	return result;
+}
 
 static std::random_device rd;
 static std::mt19937 rng(rd());
-static std::uniform_real_distribution<float> uniformDist;
 static std::normal_distribution<float> normalDist;
 
+
+void randomSeedTime() {
+	struct timeval tv;
+	gettimeofday(&tv, NULL);
+	xoroshiro128plus_state[0] = tv.tv_sec;
+	xoroshiro128plus_state[1] = tv.tv_usec;
+	// Generate a few times to fix the fact that the time is not a uniform u64
+	for (int i = 0; i < 10; i++) {
+		xoroshiro128plus_next();
+	}
+}
+
 uint32_t randomu32() {
-	return rng();
+	return xoroshiro128plus_next() >> 32;
+}
+
+uint64_t randomu64() {
+	return xoroshiro128plus_next();
 }
 
 float randomf() {
-	return uniformDist(rng);
+	// 24 bits of granularity is the best that can be done with floats while ensuring that the return value lies in [0.0, 1.0).
+	return (xoroshiro128plus_next() >> (64 - 24)) / powf(2, 24);
 }
 
-float randomNormal(){
-	return normalDist(rng);
+float randomNormal() {
+	// Box-Muller transform
+	float radius = sqrtf(-2.f * logf(1.f - randomf()));
+	float theta = 2.f * M_PI * randomf();
+	return radius * sinf(theta);
+
+	// // Central Limit Theorem
+	// const int n = 8;
+	// float sum = 0.0;
+	// for (int i = 0; i < n; i++) {
+	// 	sum += randomf();
+	// }
+	// return (sum - n / 2.f) / sqrtf(n / 12.f);
 }
 
 
@@ -45,7 +92,7 @@ std::string stringf(const char *format, ...) {
 	std::string s;
 	s.resize(size);
 	va_start(args, format);
-	vsnprintf(&s[0], size+1, format, args);
+	vsnprintf(&s[0], size + 1, format, args);
 	va_end(args);
 	return s;
 }