Re: [PATCH] PR fortran/91414: Improved PRNG

2019-08-13 Thread Janne Blomqvist
On Mon, Aug 12, 2019 at 11:23 PM Steve Kargl
 wrote:
>
> On Sun, Aug 11, 2019 at 12:37:49PM +0300, Janne Blomqvist wrote:
> > Update the PRNG from xorshift1024* to xoshiro256** by the same
> > author. For details see
> >
> > http://prng.di.unimi.it/
> >
> > and the paper at
> >
> > https://arxiv.org/abs/1805.01407
> >
> > Also the seeding is slightly improved, by reading only 8 bytes from
> > the operating system and using the simple splitmix64 PRNG to fill in
> > the rest of the PRNG state (as recommended by the xoshiro author),
> > instead of reading the entire state from the OS.
> >
> > Regtested on x86_64-pc-linux-gnu, Ok for trunk?
> >
>
> Looks good to me.

Thanks, committed to trunk (with a minor bugfix I noticed). I
backported the initialization changes to the gcc9 and gcc8 branches as
well with the attached patch.



-- 
Janne Blomqvist
From 16abae420bcff75950868a30327bf7b242e0388e Mon Sep 17 00:00:00 2001
From: Janne Blomqvist 
Date: Tue, 13 Aug 2019 10:35:05 +0300
Subject: [PATCH] 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  

	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.
---
 libgfortran/intrinsics/random.c | 51 ++---
 1 file changed, 21 insertions(+), 30 deletions(-)

diff --git a/libgfortran/intrinsics/random.c b/libgfortran/intrinsics/random.c
index 7476439647c..9283477482f 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,13 @@ 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++)
+	{
+  os_seed = splitmix64 (os_seed);
+  master_state[i] = os_seed;
+}
   njumps = 0;
   master_init = true;
 }
-- 
2.17.1



Re: [PATCH] PR fortran/91414: Improved PRNG

2019-08-12 Thread Steve Kargl
On Sun, Aug 11, 2019 at 12:37:49PM +0300, Janne Blomqvist wrote:
> Update the PRNG from xorshift1024* to xoshiro256** by the same
> author. For details see
> 
> http://prng.di.unimi.it/
> 
> and the paper at
> 
> https://arxiv.org/abs/1805.01407
> 
> Also the seeding is slightly improved, by reading only 8 bytes from
> the operating system and using the simple splitmix64 PRNG to fill in
> the rest of the PRNG state (as recommended by the xoshiro author),
> instead of reading the entire state from the OS.
> 
> Regtested on x86_64-pc-linux-gnu, Ok for trunk?
> 

Looks good to me.

-- 
Steve


[PATCH] PR fortran/91414: Improved PRNG

2019-08-11 Thread Janne Blomqvist
Update the PRNG from xorshift1024* to xoshiro256** by the same
author. For details see

http://prng.di.unimi.it/

and the paper at

https://arxiv.org/abs/1805.01407

Also the seeding is slightly improved, by reading only 8 bytes from
the operating system and using the simple splitmix64 PRNG to fill in
the rest of the PRNG state (as recommended by the xoshiro author),
instead of reading the entire state from the OS.

Regtested on x86_64-pc-linux-gnu, Ok for trunk?

gcc/fortran/ChangeLog:

2019-08-11  Janne Blomqvist  

PR fortran/91414
* check.c (gfc_check_random_seed): Reduce seed_size.
* intrinsic.texi (RANDOM_NUMBER): Update to match new PRNG.

gcc/testsuite/ChangeLog:

2019-08-11  Janne Blomqvist  

PR fortran/91414
* gfortran.dg/random_seed_1.f90: Update to match new seed size.

libgfortran/ChangeLog:

2019-08-11  Janne Blomqvist  

PR fortran/91414
* intrinsics/random.c (prng_state): Update state struct.
(master_state): Update to match new size.
(get_rand_state): Update to match new PRNG.
(rotl): New function.
(xorshift1024star): Replace with prng_next.
(prng_next): New function.
(jump): Update for new PRNG.
(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.
(random_r4): Update to new function and struct names.
(random_r8): Likewise.
(random_r10): Likewise.
(random_r16): Likewise.
(arandom_r4): Liekwise.
(arandom_r8): Likewise.
(arandom_r10): Likwewise.
(arandom_r16): Likewise.
(xor_keys): Reduce size to match new PRNG.
(random_seed_i4): Update to new function and struct names, remove
special handling of variable p used in previous PRNG.
(random_seed_i8): Likewise.
---
 gcc/fortran/check.c |   5 +-
 gcc/fortran/intrinsic.texi  |  10 +-
 gcc/testsuite/gfortran.dg/random_seed_1.f90 |   7 +-
 libgfortran/intrinsics/random.c | 213 +---
 4 files changed, 110 insertions(+), 125 deletions(-)

diff --git a/gcc/fortran/check.c b/gcc/fortran/check.c
index 370a3c819f9..2bd8bc37556 100644
--- a/gcc/fortran/check.c
+++ b/gcc/fortran/check.c
@@ -6484,9 +6484,8 @@ gfc_check_random_seed (gfc_expr *size, gfc_expr *put, 
gfc_expr *get)
   mpz_t put_size, get_size;
 
   /* Keep the number of bytes in sync with master_state in
- libgfortran/intrinsics/random.c. +1 due to the integer p which is
- part of the state too.  */
-  seed_size = 128 / gfc_default_integer_kind + 1;
+ libgfortran/intrinsics/random.c.  */
+  seed_size = 32 / gfc_default_integer_kind;
 
   if (size != NULL)
 {
diff --git a/gcc/fortran/intrinsic.texi b/gcc/fortran/intrinsic.texi
index f390761dc3d..3aa068dba9d 100644
--- a/gcc/fortran/intrinsic.texi
+++ b/gcc/fortran/intrinsic.texi
@@ -11792,10 +11792,10 @@ end program test_random_seed
 Returns a single pseudorandom number or an array of pseudorandom numbers
 from the uniform distribution over the range @math{ 0 \leq x < 1}.
 
-The runtime-library implements the xorshift1024* random number
-generator (RNG). This generator has a period of @math{2^{1024} - 1},
-and when using multiple threads up to @math{2^{512}} threads can each
-generate @math{2^{512}} random numbers before any aliasing occurs.
+The runtime-library implements the xoshiro256** pseudorandom number
+generator (PRNG). This generator has a period of @math{2^{256} - 1},
+and when using multiple threads up to @math{2^{128}} threads can each
+generate @math{2^{128}} random numbers before any aliasing occurs.
 
 Note that in a multi-threaded program (e.g. using OpenMP directives),
 each thread will have its own random number state. For details of the
@@ -11852,7 +11852,7 @@ called either without arguments or with the @var{PUT} 
argument, the
 given seed is copied into a master seed as well as the seed of the
 current thread. When a new thread uses @code{RANDOM_NUMBER} for the
 first time, the seed is copied from the master seed, and forwarded
-@math{N * 2^{512}} steps to guarantee that the random stream does not
+@math{N * 2^{128}} steps to guarantee that the random stream does not
 alias any other stream in the system, where @var{N} is the number of
 threads that have used @code{RANDOM_NUMBER} so far during the program
 execution.
diff --git a/gcc/testsuite/gfortran.dg/random_seed_1.f90 
b/gcc/testsuite/gfortran.dg/random_seed_1.f90
index 39c81ce51b7..a97e53059f8 100644
--- a/gcc/testsuite/gfortran.dg/random_seed_1.f90
+++ b/gcc/testsuite/gfortran.dg/random_seed_1.f90
@@ -10,11 +10,12 @@
 PROGRAM random_seed_1
   IMPLICIT NONE
 
-  INTEGER, PARAMETER :: nbytes = 128
+  ! Should match sizeof(master_state) in
+  ! libgfortran/intrinsics/random.c
+  INTEGER, PARAMETER :: nbytes