Move from JKISS PRNG to Xorshift128+

random.c

@@ -1,34 +1,48 @@
 #include <stdio.h>
-#include <sys/types.h>
-#include <sys/time.h>
 #include <fcntl.h>
 #include <stdlib.h>
-#include <sys/time.h>
-#include <limits.h>
 #include <unistd.h>
 
 #include "random.h"
 #include "config.h"
 
 /* This file handles everything concerning the PRNG,
- * Currently it is JKISS by David Jones */
-
-/* Public domain code for JKISS RNG */ 
-static unsigned int x = 123456789, y = 987654321, z = 43219876, c = 6543217; /* Seed variables */ 
-static unsigned int x_save, y_save, z_save, c_save;
+ * Currently it is Xorshift128+ */
+
+/* State variables for Xorshift128+ */
+static unsigned long long s[2];
+static unsigned long long s_save[2];
+
+unsigned long long xorshift128plus() {
+    unsigned long long x = s[0];
+    const unsigned long long y = s[1];
+    s[0] = y;
+    x ^= x << 23; // a
+    s[1] = x ^ y ^ (x >> 17) ^ (y >> 26); // b, c
+    return s[1] + y;
+}
 
-unsigned int JKISS() 
-{ 
-    unsigned long long t; 
+void init_random() {
+    s[0] = ((unsigned long long)secure_rand() << 32) | secure_rand();
+    s[1] = ((unsigned long long)secure_rand() << 32) | secure_rand();
+}
 
-    x = 314527869 * x + 1234567; 
-    y ^= y << 5; y ^= y >> 7; y ^= y << 22; 
-    t = 4294584393ULL * z + c; c = t >> 32; z = t; 
+/* Called when random requires a new seed, clear_random_seed must be called
+ * to restore the previous seed values.
+ */
+void set_random_seed(unsigned int seed) {
+    s_save[0] = s[0];
+    s_save[1] = s[1];
+    s[0] = 6364136223846793005ULL * (seed ^ (seed >> 30)) + 1;
+    s[1] = 1442695040888963407ULL * (seed ^ (seed >> 27)) + 1;
+}
 
-    return x + y + z; 
-} 
+void clear_random_seed() {
+    s[0] = s_save[0];
+    s[1] = s_save[1];
+}
 
-unsigned int 
+unsigned int
 secure_rand(void) 
 { 
     int fn; 
@@ -46,31 +60,7 @@ secure_rand(void)
     return r; 
 } 
 
-void 
-init_random() 
-{ 
-    x = secure_rand(); 
-    while (!(y = secure_rand())); /* y must not be zero! */ 
-    z = secure_rand(); 
-
-    /* We don’t really need to set c as well but let's anyway… */ 
-    /* NOTE: offset c by 1 to avoid z=c=0 */ 
-    c = secure_rand() % 698769068 + 1; /* Should be less than 698769069 */ 
-} 
-
-/* Called when random requires a new seed, clear_random_seed must be called
- * to restore the previous seed values.
- */
-void
-set_random_seed(unsigned int seed) {
-    x_save = x; y_save = y; z_save = z; c_save = c;
-    x = 123456789; y = seed; z = 43219876; c = 6543217;
-}
 
-void
-clear_random_seed() {
-    x = x_save; y = y_save; z = z_save; c = c_save;
-}
 /*
  * Return a random argument in the range 0 .. n-1.
  */
@@ -79,17 +69,12 @@ random_number(long long n)
 {
     if (n <= 0)
         return 0;
-    return ((unsigned long long)JKISS() * 4294967296 + (unsigned long long)JKISS()) % n;
+    return xorshift128plus() % n;
 }
 
 double
 random_double()
 {
-    double x; 
-    unsigned int a, b; 
-
-    a = JKISS() >> 6; /* Upper 26 bits */ 
-    b = JKISS() >> 5; /* Upper 27 bits */ 
-    x = (a * 134217728.0 + b) / 9007199254740992.0; 
-    return x; 
+    unsigned long long r = xorshift128plus();
+    return (double)r / (double)0xFFFFFFFFFFFFFFFFULL;
 }

regress/tests/test-105-random.c

@@ -0,0 +1,41 @@
+#pragma strict_types
+
+void
+create()
+{
+    mapping res = ([]);
+    for (int i = 0; i < 100000; i++)
+    {
+        res[random(100000)]++;
+    }
+
+    if (m_sizeof(res) < 60000 || m_sizeof(res) > 70000)
+        throw("random() is not random enough " + m_sizeof(res) + "\n");
+    int tripped = 0;
+    foreach (int i, int j : res)
+    {
+        if (j > 7) {
+            tripped++;
+        }
+    }
+    if (tripped > 5)
+        throw("random() is not random enough tripped = " + tripped + "\n");
+
+    for (int i = 0; i < 100; i++)
+    {
+        int seed = random(100000000000000);
+        int first = random(10000000, seed);
+        for (int j = 0; j < 100; j++)
+            if (random(10000000, seed) != first)
+                write("random() not keep seeds\n");
+    }
+
+    for (int i = 0; i < 10000; i++)
+    {
+        float rand = rnd();
+        if (rand < 0.0 || rand > 1.0)
+        {
+            throw("rnd() not in range 0.0-1.0\n");
+        }
+    }
+}