Replace KISS PRNG with xorshift1024* using per-thread state.

frontend:

2016-08-11  Janne Blomqvist  <jb@gcc.gnu.org>

	* check.c (gfc_check_random_seed): Use new seed size in check.
	* intrinsic.texi (RANDOM_NUMBER): Updated documentation.
	(RANDOM_SEED): Likewise.


testsuite:

2016-08-11  Janne Blomqvist  <jb@gcc.gnu.org>

	* gfortran.dg/random_7.f90: Take into account that the last seed
	value is the special p value.
	* gfortran.dg/random_seed_1.f90: Seed size is now constant.


libgfortran:
2016-08-11  Janne Blomqvist  <jb@gcc.gnu.org>

	* intrinsics/random.c: Replace KISS with xorshift1024* using
	per-thread state.
	* runtime/main.c (init): Don't call random_seed_i4.

From-SVN: r239356

6 files changed, 57 insertions, 96 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index ef8fc17..4b9b739 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,3 +1,9 @@
+2016-08-11  Janne Blomqvist  <jb@gcc.gnu.org>
+
+	* check.c (gfc_check_random_seed): Use new seed size in check.
+	* intrinsic.texi (RANDOM_NUMBER): Updated documentation.
+	(RANDOM_SEED): Likewise.
+
 2016-08-08  Jakub Jelinek  <jakub@redhat.com>
 
 	PR fortran/72716
diff --git a/gcc/fortran/check.c b/gcc/fortran/check.c
index 288957a..ff5e80b 100644
--- a/gcc/fortran/check.c
+++ b/gcc/fortran/check.c
@@ -5527,16 +5527,14 @@ gfc_check_random_number (gfc_expr *harvest)
 bool
 gfc_check_random_seed (gfc_expr *size, gfc_expr *put, gfc_expr *get)
 {
-  unsigned int nargs = 0, kiss_size;
+  unsigned int nargs = 0, seed_size;
   locus *where = NULL;
   mpz_t put_size, get_size;
-  bool have_gfc_real_16; /* Try and mimic HAVE_GFC_REAL_16 in libgfortran.  */
 
-  have_gfc_real_16 = gfc_validate_kind (BT_REAL, 16, true) != -1;
-
-  /* Keep the number of bytes in sync with kiss_size in
-     libgfortran/intrinsics/random.c.  */
-  kiss_size = (have_gfc_real_16 ? 48 : 32) / gfc_default_integer_kind;
+  /* 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;
 
   if (size != NULL)
     {
@@ -5579,11 +5577,11 @@ gfc_check_random_seed (gfc_expr *size, gfc_expr *put, gfc_expr *get)
 	return false;
 
       if (gfc_array_size (put, &put_size)
-	  && mpz_get_ui (put_size) < kiss_size)
+	  && mpz_get_ui (put_size) < seed_size)
 	gfc_error ("Size of %qs argument of %qs intrinsic at %L "
 		   "too small (%i/%i)",
 		   gfc_current_intrinsic_arg[1]->name, gfc_current_intrinsic,
-		   where, (int) mpz_get_ui (put_size), kiss_size);
+		   where, (int) mpz_get_ui (put_size), seed_size);
     }
 
   if (get != NULL)
@@ -5611,11 +5609,11 @@ gfc_check_random_seed (gfc_expr *size, gfc_expr *put, gfc_expr *get)
 	return false;
 
        if (gfc_array_size (get, &get_size)
- 	  && mpz_get_ui (get_size) < kiss_size)
+	   && mpz_get_ui (get_size) < seed_size)
 	gfc_error ("Size of %qs argument of %qs intrinsic at %L "
 		   "too small (%i/%i)",
 		   gfc_current_intrinsic_arg[2]->name, gfc_current_intrinsic,
-		   where, (int) mpz_get_ui (get_size), kiss_size);
+		   where, (int) mpz_get_ui (get_size), seed_size);
     }
 
   /* RANDOM_SEED may not have more than one non-optional argument.  */
diff --git a/gcc/fortran/intrinsic.texi b/gcc/fortran/intrinsic.texi
index ae5d814..67b0231 100644
--- a/gcc/fortran/intrinsic.texi
+++ b/gcc/fortran/intrinsic.texi
@@ -11126,23 +11126,16 @@ end program test_rand
 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 George Marsaglia's KISS (Keep It Simple 
-Stupid) random number generator (RNG). This RNG combines:
-@enumerate
-@item The congruential generator @math{x(n) = 69069 \cdot x(n-1) + 1327217885}
-with a period of @math{2^{32}},
-@item A 3-shift shift-register generator with a period of @math{2^{32} - 1},
-@item  Two 16-bit multiply-with-carry generators with a period of
-@math{597273182964842497 > 2^{59}}.
-@end enumerate
-The overall period exceeds @math{2^{123}}.
+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.
+
+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
+seeding procedure, see the documentation for the @code{RANDOM_SEED}
+intrinsic.
 
-Please note, this RNG is thread safe if used within OpenMP directives,
-i.e., its state will be consistent while called from multiple threads.
-However, the KISS generator does not create random numbers in parallel 
-from multiple sources, but in sequence from a single source. If an
-OpenMP-enabled application heavily relies on random numbers, one should 
-consider employing a dedicated parallel random number generator instead.
 
 @item @emph{Standard}:
 Fortran 95 and later
@@ -11184,12 +11177,23 @@ end program
 Restarts or queries the state of the pseudorandom number generator used by 
 @code{RANDOM_NUMBER}.
 
-If @code{RANDOM_SEED} is called without arguments, it is initialized
-to a default state. The example below shows how to initialize the
-random seed with a varying seed in order to ensure a different random
-number sequence for each invocation of the program. Note that setting
-any of the seed values to zero should be avoided as it can result in
-poor quality random numbers being generated.
+If @code{RANDOM_SEED} is called without arguments, it is seeded with
+random data retrieved from the operating system.
+
+As an extension to the Fortran standard, the GFortran
+@code{RANDOM_NUMBER} supports multiple threads. Each thread in a
+multi-threaded program has its own seed.  When @code{RANDOM_SEED} is
+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
+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.
+
+Note that setting any of the seed values to zero should be avoided as
+it can result in poor quality random numbers being generated.
 
 @item @emph{Standard}:
 Fortran 95 and later
@@ -11217,57 +11221,16 @@ the @var{SIZE} argument.
 
 @item @emph{Example}:
 @smallexample
-subroutine init_random_seed()
-  use iso_fortran_env, only: int64
+program test_random_seed
   implicit none
   integer, allocatable :: seed(:)
-  integer :: i, n, un, istat, dt(8), pid
-  integer(int64) :: t
+  integer :: n
 
   call random_seed(size = n)
   allocate(seed(n))
-  ! First try if the OS provides a random number generator
-  open(newunit=un, file="/dev/urandom", access="stream", &
-       form="unformatted", action="read", status="old", iostat=istat)
-  if (istat == 0) then
-     read(un) seed
-     close(un)
-  else
-     ! Fallback to XOR:ing the current time and pid. The PID is
-     ! useful in case one launches multiple instances of the same
-     ! program in parallel.
-     call system_clock(t)
-     if (t == 0) then
-        call date_and_time(values=dt)
-        t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
-             + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
-             + dt(3) * 24_int64 * 60 * 60 * 1000 &
-             + dt(5) * 60 * 60 * 1000 &
-             + dt(6) * 60 * 1000 + dt(7) * 1000 &
-             + dt(8)
-     end if
-     pid = getpid()
-     t = ieor(t, int(pid, kind(t)))
-     do i = 1, n
-        seed(i) = lcg(t)
-     end do
-  end if
-  call random_seed(put=seed)
-contains
-  ! This simple PRNG might not be good enough for real work, but is
-  ! sufficient for seeding a better PRNG.
-  function lcg(s)
-    integer :: lcg
-    integer(int64) :: s
-    if (s == 0) then
-       s = 104729
-    else
-       s = mod(s, 4294967296_int64)
-    end if
-    s = mod(s * 279470273_int64, 4294967291_int64)
-    lcg = int(mod(s, int(huge(0), int64)), kind(0))
-  end function lcg
-end subroutine init_random_seed
+  call random_seed(get=seed)
+  write (*, *) seed
+end program test_random_seed
 @end smallexample
 
 @item @emph{See also}:
diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index 6798ff8..7bca99a 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,9 @@
+2016-08-11  Janne Blomqvist  <jb@gcc.gnu.org>
+
+	* gfortran.dg/random_7.f90: Take into account that the last seed
+	value is the special p value.
+	* gfortran.dg/random_seed_1.f90: Seed size is now constant.
+
 2016-08-11  Richard Biener  <rguenther@suse.de>
 
 	* gcc.dg/tree-ssa/ssa-dom-thread-7.c: Adjust.
diff --git a/gcc/testsuite/gfortran.dg/random_7.f90 b/gcc/testsuite/gfortran.dg/random_7.f90
index 6435a34..aafe346d 100644
--- a/gcc/testsuite/gfortran.dg/random_7.f90
+++ b/gcc/testsuite/gfortran.dg/random_7.f90
@@ -10,8 +10,8 @@ program trs
   seed(:) = huge(seed) / 17
   call test_random_seed(put=seed)
   call test_random_seed(get=check)
-  print *, seed
-  print *, check
+  ! In the current implementation seed(17) is special
+  seed(17) = check(17)
   if (any (seed /= check)) call abort
 contains
   subroutine test_random_seed(size, put, get)
diff --git a/gcc/testsuite/gfortran.dg/random_seed_1.f90 b/gcc/testsuite/gfortran.dg/random_seed_1.f90
index ccbcf00..39c81ce 100644
--- a/gcc/testsuite/gfortran.dg/random_seed_1.f90
+++ b/gcc/testsuite/gfortran.dg/random_seed_1.f90
@@ -3,10 +3,6 @@
 ! Emit a diagnostic for too small PUT array at compile time
 ! See PR fortran/37159
 
-! Possible improvement:
-! Provide a separate testcase for systems that support REAL(16),
-! to test the minimum size of 12 (instead of 8).
-!
 ! Updated to check for arrays of unexpected size,
 ! this also works for -fdefault-integer-8.
 !
@@ -14,19 +10,11 @@
 PROGRAM random_seed_1
   IMPLICIT NONE
 
-  ! Find out what the's largest kind size
-  INTEGER, PARAMETER :: k1 = kind (0.d0)
-  INTEGER, PARAMETER :: &
-    k2 = max (k1, selected_real_kind (precision (0._k1) + 1))
-  INTEGER, PARAMETER :: &
-    k3 = max (k2, selected_real_kind (precision (0._k2) + 1))
-  INTEGER, PARAMETER :: &
-    k4 = max (k3, selected_real_kind (precision (0._k3) + 1))
-
-  INTEGER, PARAMETER :: nbytes = MERGE(48, 32, k4 == 16)
+  INTEGER, PARAMETER :: nbytes = 128
 
+  ! +1 due to the special 'p' value in xorshift1024*
   ! '+1' to avoid out-of-bounds warnings
-  INTEGER, PARAMETER    :: n = nbytes / KIND(n) + 1
+  INTEGER, PARAMETER    :: n = nbytes / KIND(n) + 2
   INTEGER, DIMENSION(n) :: seed
 
   ! Get seed, array too small