PR fortran/91414 Improve initialization of PRNG

As part of PR 91414 an improved PRNG was contributed to trunk. This is
a partial backport of some related changes to the PRNG. Namely when
seeding the PRNG, it needs only 8 bytes of randomness from the OS, and
uses a simple splitmix64 PRNG to fill in the rest of the state,
instead of getting all the state from the OS. This can be useful for
operating systems that can run out of entropy.

libgfortran/ChangeLog:

2019-08-13  Janne Blomqvist  <jb@gcc.gnu.org>

	Partial backport from trunk
	PR fortran/91414
	* intrinsics/random.c (lcg_parkmiller): Replace with splitmix64.
	(splitmix64): New function.
	(getosrandom): Fix return value, simplify.
	(init_rand_state): Use getosrandom only to get 8 bytes, splitmix64
	to fill rest of state.


git-svn-id: https://gcc.gnu.org/svn/gcc/branches/gcc-9-branch@274362 138bc75d-0d04-0410-961f-82ee72b054a4

2 files changed, 28 insertions, 30 deletions

diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog
index 48f9886896e..70ed7487eb7 100644
--- a/libgfortran/ChangeLog
+++ b/libgfortran/ChangeLog
@@ -1,3 +1,13 @@
+2019-08-13  Janne Blomqvist  <jb@gcc.gnu.org>
+
+	Partial backport from trunk
+	PR fortran/91414
+	* intrinsics/random.c (lcg_parkmiller): Replace with splitmix64.
+	(splitmix64): New function.
+	(getosrandom): Fix return value, simplify.
+	(init_rand_state): Use getosrandom only to get 8 bytes, splitmix64
+	to fill rest of state.
+
 2019-08-12  Release Manager
 
 	* GCC 9.2.0 released.
diff --git a/libgfortran/intrinsics/random.c b/libgfortran/intrinsics/random.c
index 7476439647c..75d1315e450 100644
--- a/libgfortran/intrinsics/random.c
+++ b/libgfortran/intrinsics/random.c
@@ -275,30 +275,19 @@ jump (xorshift1024star_state* rs)
 }
 
 
-/* Super-simple LCG generator used in getosrandom () if /dev/urandom
-   doesn't exist.  */
+/* Splitmix64 recommended by xorshift author for initializing.  After
+   getting one uint64_t value from the OS, this is used to fill in the
+   rest of the state.  */
 
-#define M 2147483647 /* 2^31 - 1 (A large prime number) */
-#define A 16807      /* Prime root of M, passes statistical tests and produces a full cycle */
-#define Q 127773 /* M / A (To avoid overflow on A * seed) */
-#define R 2836   /* M % A (To avoid overflow on A * seed) */
-
-__attribute__((unused)) static uint32_t
-lcg_parkmiller(uint32_t seed)
+static uint64_t
+splitmix64 (uint64_t x)
 {
-    uint32_t hi = seed / Q;
-    uint32_t lo = seed % Q;
-    int32_t test = A * lo - R * hi;
-    if (test <= 0)
-        test += M;
-    return test;
+  uint64_t z = (x += 0x9e3779b97f4a7c15);
+  z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+  z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+  return z ^ (z >> 31);
 }
 
-#undef M
-#undef A
-#undef Q
-#undef R
-
 
 /* Get some random bytes from the operating system in order to seed
    the PRNG.  */
@@ -315,7 +304,7 @@ getosrandom (void *buf, size_t buflen)
 #else
 #ifdef HAVE_GETENTROPY
   if (getentropy (buf, buflen) == 0)
-    return 0;
+    return buflen;
 #endif
   int flags = O_RDONLY;
 #ifdef O_CLOEXEC
@@ -328,7 +317,7 @@ getosrandom (void *buf, size_t buflen)
       close (fd);
       return res;
     }
-  uint32_t seed = 1234567890;
+  uint64_t seed = 0x047f7684e9fc949dULL;
   time_t secs;
   long usecs;
   if (gf_gettime (&secs, &usecs) == 0)
@@ -340,13 +329,9 @@ getosrandom (void *buf, size_t buflen)
   pid_t pid = getpid();
   seed ^= pid;
 #endif
-  uint32_t* ub = buf;
-  for (size_t i = 0; i < buflen / sizeof (uint32_t); i++)
-    {
-      ub[i] = seed;
-      seed = lcg_parkmiller (seed);
-    }
-  return buflen;
+  size_t size = buflen < sizeof (uint64_t) ? buflen : sizeof (uint64_t);
+  memcpy (buf, &seed, size);
+  return size;
 #endif /* __MINGW64_VERSION_MAJOR  */
 }
 
@@ -361,7 +346,10 @@ init_rand_state (xorshift1024star_state* rs, const bool locked)
     __gthread_mutex_lock (&random_lock);
   if (!master_init)
     {
-      getosrandom (master_state, sizeof (master_state));
+      uint64_t os_seed;
+      getosrandom (&os_seed, sizeof (os_seed));
+      for (uint64_t i = 0; i < sizeof (master_state) / sizeof (uint64_t); i++)
+	master_state[i] = splitmix64 (os_seed);
       njumps = 0;
       master_init = true;
     }