From b152f5a2b33b251ab1874a43d97ce73d11eec0a4 Mon Sep 17 00:00:00 2001 From: Janne Blomqvist Date: Thu, 11 Aug 2016 11:58:55 +0300 Subject: [PATCH] Replace KISS PRNG with xorshift1024* using per-thread state. frontend: 2016-08-11 Janne Blomqvist * 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 * 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 * intrinsics/random.c: Replace KISS with xorshift1024* using per-thread state. * runtime/main.c (init): Don't call random_seed_i4. From-SVN: r239356 --- gcc/fortran/ChangeLog | 6 + gcc/fortran/check.c | 20 +- gcc/fortran/intrinsic.texi | 99 +--- gcc/testsuite/ChangeLog | 6 + gcc/testsuite/gfortran.dg/random_7.f90 | 4 +- gcc/testsuite/gfortran.dg/random_seed_1.f90 | 18 +- libgfortran/ChangeLog | 6 + libgfortran/intrinsics/random.c | 610 ++++++++++++-------- libgfortran/runtime/main.c | 2 - 9 files changed, 434 insertions(+), 337 deletions(-) diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog index ef8fc17603c..4b9b7398e80 100644 --- a/gcc/fortran/ChangeLog +++ b/gcc/fortran/ChangeLog @@ -1,3 +1,9 @@ +2016-08-11 Janne Blomqvist + + * 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 PR fortran/72716 diff --git a/gcc/fortran/check.c b/gcc/fortran/check.c index 288957aaa34..ff5e80b9df5 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 ae5d814004f..67b02315344 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 6798ff8ae48..7bca99a8b3b 100644 --- a/gcc/testsuite/ChangeLog +++ b/gcc/testsuite/ChangeLog @@ -1,3 +1,9 @@ +2016-08-11 Janne Blomqvist + + * 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 * 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 6435a34cf58..aafe346de7b 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 ccbcf00cf12..39c81ce51b7 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 diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog index 8b21527d0a4..ca2c8edc45c 100644 --- a/libgfortran/ChangeLog +++ b/libgfortran/ChangeLog @@ -1,3 +1,9 @@ +2016-08-11 Janne Blomqvist + + * intrinsics/random.c: Replace KISS with xorshift1024* using + per-thread state. + * runtime/main.c (init): Don't call random_seed_i4. + 2016-07-22 Andre Vehreschild * caf/libcaf.h: Add parameter stat to caf_get() and diff --git a/libgfortran/intrinsics/random.c b/libgfortran/intrinsics/random.c index 9b3473310c8..6084ebe4bb2 100644 --- a/libgfortran/intrinsics/random.c +++ b/libgfortran/intrinsics/random.c @@ -1,7 +1,7 @@ /* Implementation of the RANDOM intrinsics Copyright (C) 2002-2016 Free Software Foundation, Inc. - Contributed by Lars Segerlund - and Steve Kargl. + Contributed by Lars Segerlund , + Steve Kargl and Janne Blomqvist. This file is part of the GNU Fortran runtime library (libgfortran). @@ -28,6 +28,22 @@ see the files COPYING3 and COPYING.RUNTIME respectively. If not, see #include #include +/* For getosrandom. */ +#ifdef HAVE_SYS_TYPES_H +#include +#endif +#ifdef HAVE_UNISTD_H +#include +#endif +#include +#include +#include "time_1.h" + +#ifdef __MINGW32__ +#define HAVE_GETPID 1 +#include +#endif + extern void random_r4 (GFC_REAL_4 *); iexport_proto(random_r4); @@ -141,154 +157,222 @@ rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2) + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(0x1.p-128); } #endif -/* libgfortran previously had a Mersenne Twister, taken from the paper: - - Mersenne Twister: 623-dimensionally equidistributed - uniform pseudorandom generator. - - by Makoto Matsumoto & Takuji Nishimura - which appeared in the: ACM Transactions on Modelling and Computer - Simulations: Special Issue on Uniform Random Number - Generation. ( Early in 1998 ). - - The Mersenne Twister code was replaced due to - - (1) Simple user specified seeds lead to really bad sequences for - nearly 100000 random numbers. - (2) open(), read(), and close() were not properly declared via header - files. - (3) The global index i was abused and caused unexpected behavior with - GET and PUT. - (4) See PR 15619. - - - libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid) - random number generator. This PRNG combines: - - (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period - of 2^32, - (2) A 3-shift shift-register generator with a period of 2^32-1, - (3) Two 16-bit multiply-with-carry generators with a period of - 597273182964842497 > 2^59. - - The overall period exceeds 2^123. - - http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu - - The above web site has an archive of a newsgroup posting from George - Marsaglia with the statement: - - Subject: Random numbers for C: Improvements. - Date: Fri, 15 Jan 1999 11:41:47 -0500 - From: George Marsaglia - Message-ID: <369F6FCA.74C7C041@stat.fsu.edu> - References: <369B5E30.65A55FD1@stat.fsu.edu> - Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis - Lines: 93 - - As I hoped, several suggestions have led to - improvements in the code for RNG's I proposed for - use in C. (See the thread "Random numbers for C: Some - suggestions" in previous postings.) The improved code - is listed below. - - A question of copyright has also been raised. Unlike - DIEHARD, there is no copyright on the code below. You - are free to use it in any way you want, but you may - wish to acknowledge the source, as a courtesy. - -"There is no copyright on the code below." included the original -KISS algorithm. */ - -/* We use three KISS random number generators, with different - seeds. - As a matter of Quality of Implementation, the random numbers - we generate for different REAL kinds, starting from the same - seed, are always the same up to the precision of these types. - We do this by using three generators with different seeds, the - first one always for the most significant bits, the second one - for bits 33..64 (if present in the REAL kind), and the third one - (called twice) for REAL(16). */ - -#define GFC_SL(k, n) ((k)^((k)<<(n))) -#define GFC_SR(k, n) ((k)^((k)>>(n))) - -/* Reference for the seed: - From: "George Marsaglia" - Newsgroups: sci.math - Message-ID: - - The KISS RNG uses four seeds, x, y, z, c, - with 0<=x<2^32, 0p; + const uint64_t s0 = rs->s[p]; + uint64_t s1 = rs->s[p = (p + 1) & 15]; + s1 ^= s1 << 31; + rs->s[p] = s1 ^ s0 ^ (s1 >> 11) ^ (s0 >> 30); + rs->p = p; + return rs->s[p] * UINT64_C(1181783497276652981); +} - seed[0] = 69069 * seed[0] + 1327217885; - seed[1] = GFC_SL(GFC_SR(GFC_SL(seed[1],13),17),5); - seed[2] = 18000 * (seed[2] & 65535) + (seed[2] >> 16); - seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16); - kiss = seed[0] + seed[1] + (seed[2] << 16) + seed[3]; - return kiss; +/* This is the jump function for the generator. It is equivalent to + 2^512 calls to xorshift1024star(); it can be used to generate 2^512 + non-overlapping subsequences for parallel computations. */ + +static void +jump (xorshift1024star_state* rs) +{ + static const uint64_t JUMP[] = { + 0x84242f96eca9c41dULL, 0xa3c65b8776f96855ULL, 0x5b34a39f070b5837ULL, + 0x4489affce4f31a1eULL, 0x2ffeeb0a48316f40ULL, 0xdc2d9891fe68c022ULL, + 0x3659132bb12fea70ULL, 0xaac17d8efa43cab8ULL, 0xc4cb815590989b13ULL, + 0x5ee975283d71c93bULL, 0x691548c86c1bd540ULL, 0x7910c41d10a1e6a5ULL, + 0x0b5fc64563b3e2a8ULL, 0x047f7684e9fc949dULL, 0xb99181f2d8f685caULL, + 0x284600e3f30e38c3ULL + }; + + uint64_t t[16] = { 0 }; + for(unsigned int i = 0; i < sizeof JUMP / sizeof *JUMP; i++) + for(int b = 0; b < 64; b++) + { + if (JUMP[i] & 1ULL << b) + for(int j = 0; j < 16; j++) + t[j] ^= rs->s[(j + rs->p) & 15]; + xorshift1024star (rs); + } + for(int j = 0; j < 16; j++) + rs->s[(j + rs->p) & 15] = t[j]; } + +/* Initialize the random number generator for the current thread, + using the master state and the number of times we must jump. */ + +static void +init_rand_state (xorshift1024star_state* rs, const bool locked) +{ + if (!locked) + __gthread_mutex_lock (&random_lock); + memcpy (&rs->s, master_state, sizeof (master_state)); + unsigned n = njumps++; + if (!locked) + __gthread_mutex_unlock (&random_lock); + for (unsigned i = 0; i < n; i++) + jump (rs); + rs->init = true; +} + + +/* Super-simple LCG generator used in getosrandom () if /dev/urandom + doesn't exist. */ + +#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) */ + +static uint32_t +lcg_parkmiller(uint32_t seed) +{ + uint32_t hi = seed / Q; + uint32_t lo = seed % Q; + int32_t test = A * lo - R * hi; + if (test <= 0) + test += M; + return test; +} + +#undef M +#undef A +#undef Q +#undef R + +/* Get some random bytes from the operating system in order to seed + the PRNG. */ + +static int +getosrandom (void *buf, size_t buflen) +{ + /* TODO: On Windows one could use CryptGenRandom + + TODO: When glibc adds a wrapper for the getrandom() system call + on Linux, one could use that. + + TODO: One could use getentropy() on OpenBSD. */ + int flags = O_RDONLY; +#ifdef O_CLOEXEC + flags |= O_CLOEXEC; +#endif + int fd = open("/dev/urandom", flags); + if (fd != -1) + { + int res = read(fd, buf, buflen); + close (fd); + return res; + } + uint32_t seed = 1234567890; + time_t secs; + long usecs; + if (gf_gettime (&secs, &usecs) == 0) + { + seed ^= secs; + seed ^= usecs; + } +#ifdef HAVE_GETPID + pid_t pid = getpid(); + seed ^= pid; +#endif + uint32_t* ub = buf; + for (size_t i = 0; i < buflen; i++) + { + ub[i] = seed; + seed = lcg_parkmiller (seed); + } + return buflen; +} + + /* This function produces a REAL(4) value from the uniform distribution with range [0,1). */ void random_r4 (GFC_REAL_4 *x) { - GFC_UINTEGER_4 kiss; - - __gthread_mutex_lock (&random_lock); - kiss = kiss_random_kernel (kiss_seed_1); - rnumber_4 (x, kiss); - __gthread_mutex_unlock (&random_lock); + xorshift1024star_state* rs = get_rand_state(); + + if (unlikely (!rs->init)) + init_rand_state (rs, false); + uint64_t r = xorshift1024star (rs); + /* Take the higher bits, ensuring that a stream of real(4), real(8), + and real(10) will be identical (except for precision). */ + uint32_t high = (uint32_t) (r >> 32); + rnumber_4 (x, high); } iexport(random_r4); @@ -298,13 +382,13 @@ iexport(random_r4); void random_r8 (GFC_REAL_8 *x) { - GFC_UINTEGER_8 kiss; + GFC_UINTEGER_8 r; + xorshift1024star_state* rs = get_rand_state(); - __gthread_mutex_lock (&random_lock); - kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; - kiss += kiss_random_kernel (kiss_seed_2); - rnumber_8 (x, kiss); - __gthread_mutex_unlock (&random_lock); + if (unlikely (!rs->init)) + init_rand_state (rs, false); + r = xorshift1024star (rs); + rnumber_8 (x, r); } iexport(random_r8); @@ -316,13 +400,13 @@ iexport(random_r8); void random_r10 (GFC_REAL_10 *x) { - GFC_UINTEGER_8 kiss; + GFC_UINTEGER_8 r; + xorshift1024star_state* rs = get_rand_state(); - __gthread_mutex_lock (&random_lock); - kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; - kiss += kiss_random_kernel (kiss_seed_2); - rnumber_10 (x, kiss); - __gthread_mutex_unlock (&random_lock); + if (unlikely (!rs->init)) + init_rand_state (rs, false); + r = xorshift1024star (rs); + rnumber_10 (x, r); } iexport(random_r10); @@ -336,22 +420,20 @@ iexport(random_r10); void random_r16 (GFC_REAL_16 *x) { - GFC_UINTEGER_8 kiss1, kiss2; - - __gthread_mutex_lock (&random_lock); - kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; - kiss1 += kiss_random_kernel (kiss_seed_2); - - kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32; - kiss2 += kiss_random_kernel (kiss_seed_3); - - rnumber_16 (x, kiss1, kiss2); - __gthread_mutex_unlock (&random_lock); + GFC_UINTEGER_8 r1, r2; + xorshift1024star_state* rs = get_rand_state(); + + if (unlikely (!rs->init)) + init_rand_state (rs, false); + r1 = xorshift1024star (rs); + r2 = xorshift1024star (rs); + rnumber_16 (x, r1, r2); } iexport(random_r16); #endif + /* This function fills a REAL(4) array with values from the uniform distribution with range [0,1). */ @@ -364,9 +446,10 @@ arandom_r4 (gfc_array_r4 *x) index_type stride0; index_type dim; GFC_REAL_4 *dest; - GFC_UINTEGER_4 kiss; + xorshift1024star_state* rs = get_rand_state(); int n; + dest = x->base_addr; dim = GFC_DESCRIPTOR_RANK (x); @@ -382,13 +465,15 @@ arandom_r4 (gfc_array_r4 *x) stride0 = stride[0]; - __gthread_mutex_lock (&random_lock); + if (unlikely (!rs->init)) + init_rand_state (rs, false); while (dest) { /* random_r4 (dest); */ - kiss = kiss_random_kernel (kiss_seed_1); - rnumber_4 (dest, kiss); + uint64_t r = xorshift1024star (rs); + uint32_t high = (uint32_t) (r >> 32); + rnumber_4 (dest, high); /* Advance to the next element. */ dest += stride0; @@ -416,7 +501,6 @@ arandom_r4 (gfc_array_r4 *x) } } } - __gthread_mutex_unlock (&random_lock); } /* This function fills a REAL(8) array with values from the uniform @@ -431,7 +515,7 @@ arandom_r8 (gfc_array_r8 *x) index_type stride0; index_type dim; GFC_REAL_8 *dest; - GFC_UINTEGER_8 kiss; + xorshift1024star_state* rs = get_rand_state(); int n; dest = x->base_addr; @@ -449,14 +533,14 @@ arandom_r8 (gfc_array_r8 *x) stride0 = stride[0]; - __gthread_mutex_lock (&random_lock); + if (unlikely (!rs->init)) + init_rand_state (rs, false); while (dest) { /* random_r8 (dest); */ - kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; - kiss += kiss_random_kernel (kiss_seed_2); - rnumber_8 (dest, kiss); + uint64_t r = xorshift1024star (rs); + rnumber_8 (dest, r); /* Advance to the next element. */ dest += stride0; @@ -484,7 +568,6 @@ arandom_r8 (gfc_array_r8 *x) } } } - __gthread_mutex_unlock (&random_lock); } #ifdef HAVE_GFC_REAL_10 @@ -501,7 +584,7 @@ arandom_r10 (gfc_array_r10 *x) index_type stride0; index_type dim; GFC_REAL_10 *dest; - GFC_UINTEGER_8 kiss; + xorshift1024star_state* rs = get_rand_state(); int n; dest = x->base_addr; @@ -519,14 +602,14 @@ arandom_r10 (gfc_array_r10 *x) stride0 = stride[0]; - __gthread_mutex_lock (&random_lock); + if (unlikely (!rs->init)) + init_rand_state (rs, false); while (dest) { /* random_r10 (dest); */ - kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; - kiss += kiss_random_kernel (kiss_seed_2); - rnumber_10 (dest, kiss); + uint64_t r = xorshift1024star (rs); + rnumber_10 (dest, r); /* Advance to the next element. */ dest += stride0; @@ -554,7 +637,6 @@ arandom_r10 (gfc_array_r10 *x) } } } - __gthread_mutex_unlock (&random_lock); } #endif @@ -573,7 +655,7 @@ arandom_r16 (gfc_array_r16 *x) index_type stride0; index_type dim; GFC_REAL_16 *dest; - GFC_UINTEGER_8 kiss1, kiss2; + xorshift1024star_state* rs = get_rand_state(); int n; dest = x->base_addr; @@ -591,18 +673,15 @@ arandom_r16 (gfc_array_r16 *x) stride0 = stride[0]; - __gthread_mutex_lock (&random_lock); + if (unlikely (!rs->init)) + init_rand_state (rs, false); while (dest) { /* random_r16 (dest); */ - kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; - kiss1 += kiss_random_kernel (kiss_seed_2); - - kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32; - kiss2 += kiss_random_kernel (kiss_seed_3); - - rnumber_16 (dest, kiss1, kiss2); + uint64_t r1 = xorshift1024star (rs); + uint64_t r2 = xorshift1024star (rs); + rnumber_16 (dest, r1, r2); /* Advance to the next element. */ dest += stride0; @@ -630,7 +709,6 @@ arandom_r16 (gfc_array_r16 *x) } } } - __gthread_mutex_unlock (&random_lock); } #endif @@ -665,22 +743,57 @@ unscramble_seed (unsigned char *dest, unsigned char *src, int size) void random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) { - unsigned char seed[4 * KISS_SIZE]; - - __gthread_mutex_lock (&random_lock); + unsigned char seed[sizeof (master_state)]; +#define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_4)) /* Check that we only have one argument present. */ if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) runtime_error ("RANDOM_SEED should have at most one argument present."); + if (size != NULL) + *size = SZ + 1; + + xorshift1024star_state* rs = get_rand_state(); + + /* Return the seed to GET data. */ + if (get != NULL) + { + /* If the rank of the array is not 1, abort. */ + if (GFC_DESCRIPTOR_RANK (get) != 1) + runtime_error ("Array rank of GET is not 1."); + + /* If the array is too small, abort. */ + if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1) + runtime_error ("Array size of GET is too small."); + + if (!rs->init) + init_rand_state (rs, false); + + /* Unscramble the seed. */ + unscramble_seed (seed, (unsigned char *) rs->s, sizeof seed); + + /* Then copy it back to the user variable. */ + for (size_t i = 0; i < SZ ; i++) + memcpy (&(get->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]), + seed + i * sizeof(GFC_UINTEGER_4), + sizeof(GFC_UINTEGER_4)); + + /* Finally copy the value of p after the seed. */ + get->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(get, 0)] = rs->p; + } + + else + { + __gthread_mutex_lock (&random_lock); + /* From the standard: "If no argument is present, the processor assigns a processor-dependent value to the seed." */ if (size == NULL && put == NULL && get == NULL) - for (size_t i = 0; i < KISS_SIZE; i++) - kiss_seed[i] = kiss_default_seed[i]; - - if (size != NULL) - *size = KISS_SIZE; + { + getosrandom (master_state, sizeof (master_state)); + njumps = 0; + init_rand_state (rs, true); + } if (put != NULL) { @@ -689,42 +802,29 @@ random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) runtime_error ("Array rank of PUT is not 1."); /* If the array is too small, abort. */ - if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) KISS_SIZE) + if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ + 1) runtime_error ("Array size of PUT is too small."); /* We copy the seed given by the user. */ - for (size_t i = 0; i < KISS_SIZE; i++) + for (size_t i = 0; i < SZ; i++) memcpy (seed + i * sizeof(GFC_UINTEGER_4), - &(put->base_addr[(KISS_SIZE - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]), + &(put->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]), sizeof(GFC_UINTEGER_4)); /* We put it after scrambling the bytes, to paper around users who provide seeds with quality only in the lower or upper part. */ - scramble_seed ((unsigned char *) kiss_seed, seed, 4 * KISS_SIZE); - } - - /* Return the seed to GET data. */ - if (get != NULL) - { - /* If the rank of the array is not 1, abort. */ - if (GFC_DESCRIPTOR_RANK (get) != 1) - runtime_error ("Array rank of GET is not 1."); + scramble_seed ((unsigned char *) master_state, seed, sizeof seed); + njumps = 0; + init_rand_state (rs, true); - /* If the array is too small, abort. */ - if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) KISS_SIZE) - runtime_error ("Array size of GET is too small."); + rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15; + } - /* Unscramble the seed. */ - unscramble_seed (seed, (unsigned char *) kiss_seed, 4 * KISS_SIZE); - /* Then copy it back to the user variable. */ - for (size_t i = 0; i < KISS_SIZE; i++) - memcpy (&(get->base_addr[(KISS_SIZE - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]), - seed + i * sizeof(GFC_UINTEGER_4), - sizeof(GFC_UINTEGER_4)); - } __gthread_mutex_unlock (&random_lock); + } +#undef SZ } iexport(random_seed_i4); @@ -732,63 +832,95 @@ iexport(random_seed_i4); void random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get) { - __gthread_mutex_lock (&random_lock); - /* Check that we only have one argument present. */ if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) runtime_error ("RANDOM_SEED should have at most one argument present."); - /* From the standard: "If no argument is present, the processor assigns - a processor-dependent value to the seed." */ - if (size == NULL && put == NULL && get == NULL) - for (size_t i = 0; i < KISS_SIZE; i++) - kiss_seed[i] = kiss_default_seed[i]; - +#define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_8)) if (size != NULL) - *size = KISS_SIZE / 2; + *size = SZ + 1; - if (put != NULL) + xorshift1024star_state* rs = get_rand_state(); + + /* Return the seed to GET data. */ + if (get != NULL) { /* If the rank of the array is not 1, abort. */ - if (GFC_DESCRIPTOR_RANK (put) != 1) - runtime_error ("Array rank of PUT is not 1."); + if (GFC_DESCRIPTOR_RANK (get) != 1) + runtime_error ("Array rank of GET is not 1."); /* If the array is too small, abort. */ - if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) KISS_SIZE / 2) - runtime_error ("Array size of PUT is too small."); + if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1) + runtime_error ("Array size of GET is too small."); + + if (!rs->init) + init_rand_state (rs, false); /* This code now should do correct strides. */ - for (size_t i = 0; i < KISS_SIZE / 2; i++) - memcpy (&kiss_seed[2*i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]), + for (size_t i = 0; i < SZ; i++) + memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &rs->s[i], sizeof (GFC_UINTEGER_8)); + + get->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(get, 0)] = rs->p; } - /* Return the seed to GET data. */ - if (get != NULL) + else + { + __gthread_mutex_lock (&random_lock); + + /* From the standard: "If no argument is present, the processor assigns + a processor-dependent value to the seed." */ + if (size == NULL && put == NULL && get == NULL) + { + getosrandom (master_state, sizeof (master_state)); + njumps = 0; + init_rand_state (rs, true); + } + + if (put != NULL) { /* If the rank of the array is not 1, abort. */ - if (GFC_DESCRIPTOR_RANK (get) != 1) - runtime_error ("Array rank of GET is not 1."); + if (GFC_DESCRIPTOR_RANK (put) != 1) + runtime_error ("Array rank of PUT is not 1."); /* If the array is too small, abort. */ - if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) KISS_SIZE / 2) - runtime_error ("Array size of GET is too small."); + if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ + 1) + runtime_error ("Array size of PUT is too small."); /* This code now should do correct strides. */ - for (size_t i = 0; i < KISS_SIZE / 2; i++) - memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &kiss_seed[2*i], + for (size_t i = 0; i < SZ; i++) + memcpy (&master_state[i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]), sizeof (GFC_UINTEGER_8)); - } + + njumps = 0; + init_rand_state (rs, true); + rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15; + } + __gthread_mutex_unlock (&random_lock); + } } iexport(random_seed_i8); -#ifndef __GTHREAD_MUTEX_INIT +#if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS static void __attribute__((constructor)) -init (void) +constructor_random (void) { +#ifndef __GTHREAD_MUTEX_INIT __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock); +#endif + if (__gthread_active_p ()) + __gthread_key_create (&rand_state_key, &free); +} +#endif + +#ifdef __GTHREADS +static void __attribute__((destructor)) +destructor_random (void) +{ + if (__gthread_active_p ()) + __gthread_key_delete (rand_state_key); } #endif diff --git a/libgfortran/runtime/main.c b/libgfortran/runtime/main.c index 5065d9cf9b3..09b89bc3df3 100644 --- a/libgfortran/runtime/main.c +++ b/libgfortran/runtime/main.c @@ -119,8 +119,6 @@ init (void) set_fpu (); init_compile_options (); - - random_seed_i4 (NULL, NULL, NULL); } -- 2.30.2