git
/
gcc
mirror of git://gcc.gnu.org/git/gcc.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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
This commit is contained in:
Janne Blomqvist 2016-08-11 11:58:55 +03:00
parent bb7ebad1c0
commit b152f5a2b3
9 changed files with 435 additions and 338 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
gcc/fortran/ChangeLog
View File

@ -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> 2016-08-08 Jakub Jelinek <jakub@redhat.com>
PR fortran/72716 PR fortran/72716

20
gcc/fortran/check.c
View File

@ -5527,16 +5527,14 @@ gfc_check_random_number (gfc_expr *harvest)
bool bool
gfc_check_random_seed (gfc_expr *size, gfc_expr *put, gfc_expr *get) 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; locus *where = NULL;
mpz_t put_size, get_size; 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 master_state in
libgfortran/intrinsics/random.c. +1 due to the integer p which is
/* Keep the number of bytes in sync with kiss_size in part of the state too. */
libgfortran/intrinsics/random.c. */ seed_size = 128 / gfc_default_integer_kind + 1;
kiss_size = (have_gfc_real_16 ? 48 : 32) / gfc_default_integer_kind;
if (size != NULL) if (size != NULL)
{ {
@ -5579,11 +5577,11 @@ gfc_check_random_seed (gfc_expr *size, gfc_expr *put, gfc_expr *get)
return false; return false;
if (gfc_array_size (put, &put_size) 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 " gfc_error ("Size of %qs argument of %qs intrinsic at %L "
"too small (%i/%i)", "too small (%i/%i)",
gfc_current_intrinsic_arg[1]->name, gfc_current_intrinsic, 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) if (get != NULL)
@ -5611,11 +5609,11 @@ gfc_check_random_seed (gfc_expr *size, gfc_expr *put, gfc_expr *get)
return false; return false;
if (gfc_array_size (get, &get_size) 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 " gfc_error ("Size of %qs argument of %qs intrinsic at %L "
"too small (%i/%i)", "too small (%i/%i)",
gfc_current_intrinsic_arg[2]->name, gfc_current_intrinsic, 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. */ /* RANDOM_SEED may not have more than one non-optional argument. */

99
gcc/fortran/intrinsic.texi
View File

@ -11126,23 +11126,16 @@ end program test_rand
Returns a single pseudorandom number or an array of pseudorandom numbers Returns a single pseudorandom number or an array of pseudorandom numbers
from the uniform distribution over the range @math{ 0 \leq x < 1}. from the uniform distribution over the range @math{ 0 \leq x < 1}.
The runtime-library implements George Marsaglia's KISS (Keep It Simple The runtime-library implements the xorshift1024* random number
Stupid) random number generator (RNG). This RNG combines: generator (RNG). This generator has a period of @math{2^{1024} - 1},
@enumerate and when using multiple threads up to @math{2^{512}} threads can each
@item The congruential generator @math{x(n) = 69069 \cdot x(n-1) + 1327217885} generate @math{2^{512}} random numbers before any aliasing occurs.
with a period of @math{2^{32}},
@item A 3-shift shift-register generator with a period of @math{2^{32} - 1}, Note that in a multi-threaded program (e.g. using OpenMP directives),
@item Two 16-bit multiply-with-carry generators with a period of each thread will have its own random number state. For details of the
@math{597273182964842497 > 2^{59}}. seeding procedure, see the documentation for the @code{RANDOM_SEED}
@end enumerate intrinsic.
The overall period exceeds @math{2^{123}}.
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}: @item @emph{Standard}:
Fortran 95 and later Fortran 95 and later
@ -11184,12 +11177,23 @@ end program
Restarts or queries the state of the pseudorandom number generator used by Restarts or queries the state of the pseudorandom number generator used by
@code{RANDOM_NUMBER}. @code{RANDOM_NUMBER}.
If @code{RANDOM_SEED} is called without arguments, it is initialized If @code{RANDOM_SEED} is called without arguments, it is seeded with
to a default state. The example below shows how to initialize the random data retrieved from the operating system.
random seed with a varying seed in order to ensure a different random
number sequence for each invocation of the program. Note that setting As an extension to the Fortran standard, the GFortran
any of the seed values to zero should be avoided as it can result in @code{RANDOM_NUMBER} supports multiple threads. Each thread in a
poor quality random numbers being generated. 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}: @item @emph{Standard}:
Fortran 95 and later Fortran 95 and later
@ -11217,57 +11221,16 @@ the @var{SIZE} argument.
@item @emph{Example}: @item @emph{Example}:
@smallexample @smallexample
subroutine init_random_seed() program test_random_seed
use iso_fortran_env, only: int64
implicit none implicit none
integer, allocatable :: seed(:) integer, allocatable :: seed(:)
integer :: i, n, un, istat, dt(8), pid integer :: n
integer(int64) :: t
call random_seed(size = n) call random_seed(size = n)
allocate(seed(n)) allocate(seed(n))
! First try if the OS provides a random number generator call random_seed(get=seed)
open(newunit=un, file="/dev/urandom", access="stream", & write (*, *) seed
form="unformatted", action="read", status="old", iostat=istat) end program test_random_seed
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
@end smallexample @end smallexample
@item @emph{See also}: @item @emph{See also}:

6
gcc/testsuite/ChangeLog
View File

@ -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> 2016-08-11 Richard Biener <rguenther@suse.de>
* gcc.dg/tree-ssa/ssa-dom-thread-7.c: Adjust. * gcc.dg/tree-ssa/ssa-dom-thread-7.c: Adjust.

4
gcc/testsuite/gfortran.dg/random_7.f90
View File

@ -10,8 +10,8 @@ program trs
seed(:) = huge(seed) / 17 seed(:) = huge(seed) / 17
call test_random_seed(put=seed) call test_random_seed(put=seed)
call test_random_seed(get=check) call test_random_seed(get=check)
print *, seed ! In the current implementation seed(17) is special
print *, check seed(17) = check(17)
if (any (seed /= check)) call abort if (any (seed /= check)) call abort
contains contains
subroutine test_random_seed(size, put, get) subroutine test_random_seed(size, put, get)

18
gcc/testsuite/gfortran.dg/random_seed_1.f90
View File

@ -3,10 +3,6 @@
! Emit a diagnostic for too small PUT array at compile time ! Emit a diagnostic for too small PUT array at compile time
! See PR fortran/37159 ! 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, ! Updated to check for arrays of unexpected size,
! this also works for -fdefault-integer-8. ! this also works for -fdefault-integer-8.
! !
@ -14,19 +10,11 @@
PROGRAM random_seed_1 PROGRAM random_seed_1
IMPLICIT NONE IMPLICIT NONE
! Find out what the's largest kind size INTEGER, PARAMETER :: nbytes = 128
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)
! +1 due to the special 'p' value in xorshift1024*
! '+1' to avoid out-of-bounds warnings ! '+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 INTEGER, DIMENSION(n) :: seed
! Get seed, array too small ! Get seed, array too small

6
libgfortran/ChangeLog
View File

@ -1,3 +1,9 @@
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.
2016-07-22 Andre Vehreschild <vehre@gcc.gnu.org> 2016-07-22 Andre Vehreschild <vehre@gcc.gnu.org>
* caf/libcaf.h: Add parameter stat to caf_get() and * caf/libcaf.h: Add parameter stat to caf_get() and

612
libgfortran/intrinsics/random.c
View File

@ -1,7 +1,7 @@
/* Implementation of the RANDOM intrinsics /* Implementation of the RANDOM intrinsics
Copyright (C) 2002-2016 Free Software Foundation, Inc. Copyright (C) 2002-2016 Free Software Foundation, Inc.
Contributed by Lars Segerlund <seger@linuxmail.org> Contributed by Lars Segerlund <seger@linuxmail.org>,
and Steve Kargl. Steve Kargl and Janne Blomqvist.
This file is part of the GNU Fortran runtime library (libgfortran). 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 <gthr.h> #include <gthr.h>
#include <string.h> #include <string.h>
/* For getosrandom. */
#ifdef HAVE_SYS_TYPES_H
#include <sys/types.h>
#endif
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
#include <sys/stat.h>
#include <fcntl.h>
#include "time_1.h"
#ifdef __MINGW32__
#define HAVE_GETPID 1
#include <process.h>
#endif
extern void random_r4 (GFC_REAL_4 *); extern void random_r4 (GFC_REAL_4 *);
iexport_proto(random_r4); iexport_proto(random_r4);
@ -141,141 +157,206 @@ 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); + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(0x1.p-128);
} }
#endif #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 We use the xorshift1024* generator, a fast high-quality generator
of 2^32, that:
(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. - passes TestU1 without any failures
http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu - provides a "jump" function making it easy to provide many
independent parallel streams.
The above web site has an archive of a newsgroup posting from George - Long period of 2**1024 - 1
Marsaglia with the statement:
Subject: Random numbers for C: Improvements. A description can be found at
Date: Fri, 15 Jan 1999 11:41:47 -0500
From: George Marsaglia <geo@stat.fsu.edu>
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 http://vigna.di.unimi.it/ftp/papers/xorshift.pdf
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 or
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 http://arxiv.org/abs/1402.6246
KISS algorithm. */
/* We use three KISS random number generators, with different The paper includes public domain source code which is the basis for
seeds. the implementation below.
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))) typedef struct
/* Reference for the seed:
From: "George Marsaglia" <g...@stat.fsu.edu>
Newsgroups: sci.math
Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
The KISS RNG uses four seeds, x, y, z, c,
with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
except that the two pairs
z=0,c=0 and z=2^32-1,c=698769068
should be avoided. */
/* Any modifications to the seeds that change KISS_SIZE below need to be
reflected in check.c (gfc_check_random_seed) to enable correct
compile-time checking of PUT size for the RANDOM_SEED intrinsic. */
#define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
#define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
#ifdef HAVE_GFC_REAL_16
#define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
#endif
static GFC_UINTEGER_4 kiss_seed[] = {
KISS_DEFAULT_SEED_1,
KISS_DEFAULT_SEED_2,
#ifdef HAVE_GFC_REAL_16
KISS_DEFAULT_SEED_3
#endif
};
static GFC_UINTEGER_4 kiss_default_seed[] = {
KISS_DEFAULT_SEED_1,
KISS_DEFAULT_SEED_2,
#ifdef HAVE_GFC_REAL_16
KISS_DEFAULT_SEED_3
#endif
};
#define KISS_SIZE (sizeof(kiss_seed)/sizeof(kiss_seed[0]))
static GFC_UINTEGER_4 * const kiss_seed_1 = kiss_seed;
static GFC_UINTEGER_4 * const kiss_seed_2 = kiss_seed + 4;
#ifdef HAVE_GFC_REAL_16
static GFC_UINTEGER_4 * const kiss_seed_3 = kiss_seed + 8;
#endif
/* kiss_random_kernel() returns an integer value in the range of
(0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
should be uniform. */
static GFC_UINTEGER_4
kiss_random_kernel(GFC_UINTEGER_4 * seed)
{ {
GFC_UINTEGER_4 kiss; bool init;
int p;
seed[0] = 69069 * seed[0] + 1327217885; uint64_t s[16];
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;
} }
xorshift1024star_state;
/* How many times we have jumped. This and master_state are the only
variables protected by random_lock. */
static unsigned njumps;
static uint64_t master_state[] = {
0xad63fa1ed3b55f36ULL, 0xd94473e78978b497ULL, 0xbc60592a98172477ULL,
0xa3de7c6e81265301ULL, 0x586640c5e785af27ULL, 0x7a2a3f63b67ce5eaULL,
0x9fde969f922d9b82ULL, 0xe6fe34379b3f3822ULL, 0x6c277eac3e99b6c2ULL,
0x9197290ab0d3f069ULL, 0xdb227302f6c25576ULL, 0xee0209aee527fae9ULL,
0x675666a793cd05b9ULL, 0xd048c99fbc70c20fULL, 0x775f8c3dba385ef5ULL,
0x625288bc262faf33ULL
};
static __gthread_key_t rand_state_key;
static xorshift1024star_state*
get_rand_state (void)
{
/* For single threaded apps. */
static xorshift1024star_state rand_state;
if (__gthread_active_p ())
{
void* p = __gthread_getspecific (rand_state_key);
if (!p)
{
p = xcalloc (1, sizeof (xorshift1024star_state));
__gthread_setspecific (rand_state_key, p);
}
return p;
}
else
return &rand_state;
}
static uint64_t
xorshift1024star (xorshift1024star_state* rs)
{
int p = rs->p;
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);
}
/* 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 /* This function produces a REAL(4) value from the uniform distribution
with range [0,1). */ with range [0,1). */
@ -283,12 +364,15 @@ kiss_random_kernel(GFC_UINTEGER_4 * seed)
void void
random_r4 (GFC_REAL_4 *x) random_r4 (GFC_REAL_4 *x)
{ {
GFC_UINTEGER_4 kiss; xorshift1024star_state* rs = get_rand_state();
__gthread_mutex_lock (&random_lock); if (unlikely (!rs->init))
kiss = kiss_random_kernel (kiss_seed_1); init_rand_state (rs, false);
rnumber_4 (x, kiss); uint64_t r = xorshift1024star (rs);
__gthread_mutex_unlock (&random_lock); /* 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); iexport(random_r4);
@ -298,13 +382,13 @@ iexport(random_r4);
void void
random_r8 (GFC_REAL_8 *x) 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); if (unlikely (!rs->init))
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; init_rand_state (rs, false);
kiss += kiss_random_kernel (kiss_seed_2); r = xorshift1024star (rs);
rnumber_8 (x, kiss); rnumber_8 (x, r);
__gthread_mutex_unlock (&random_lock);
} }
iexport(random_r8); iexport(random_r8);
@ -316,13 +400,13 @@ iexport(random_r8);
void void
random_r10 (GFC_REAL_10 *x) 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); if (unlikely (!rs->init))
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; init_rand_state (rs, false);
kiss += kiss_random_kernel (kiss_seed_2); r = xorshift1024star (rs);
rnumber_10 (x, kiss); rnumber_10 (x, r);
__gthread_mutex_unlock (&random_lock);
} }
iexport(random_r10); iexport(random_r10);
@ -336,22 +420,20 @@ iexport(random_r10);
void void
random_r16 (GFC_REAL_16 *x) random_r16 (GFC_REAL_16 *x)
{ {
GFC_UINTEGER_8 kiss1, kiss2; GFC_UINTEGER_8 r1, r2;
xorshift1024star_state* rs = get_rand_state();
__gthread_mutex_lock (&random_lock); if (unlikely (!rs->init))
kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; init_rand_state (rs, false);
kiss1 += kiss_random_kernel (kiss_seed_2); r1 = xorshift1024star (rs);
r2 = xorshift1024star (rs);
kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32; rnumber_16 (x, r1, r2);
kiss2 += kiss_random_kernel (kiss_seed_3);
rnumber_16 (x, kiss1, kiss2);
__gthread_mutex_unlock (&random_lock);
} }
iexport(random_r16); iexport(random_r16);
#endif #endif
/* This function fills a REAL(4) array with values from the uniform /* This function fills a REAL(4) array with values from the uniform
distribution with range [0,1). */ distribution with range [0,1). */
@ -364,9 +446,10 @@ arandom_r4 (gfc_array_r4 *x)
index_type stride0; index_type stride0;
index_type dim; index_type dim;
GFC_REAL_4 *dest; GFC_REAL_4 *dest;
GFC_UINTEGER_4 kiss; xorshift1024star_state* rs = get_rand_state();
int n; int n;
dest = x->base_addr; dest = x->base_addr;
dim = GFC_DESCRIPTOR_RANK (x); dim = GFC_DESCRIPTOR_RANK (x);
@ -382,13 +465,15 @@ arandom_r4 (gfc_array_r4 *x)
stride0 = stride[0]; stride0 = stride[0];
__gthread_mutex_lock (&random_lock); if (unlikely (!rs->init))
init_rand_state (rs, false);
while (dest) while (dest)
{ {
/* random_r4 (dest); */ /* random_r4 (dest); */
kiss = kiss_random_kernel (kiss_seed_1); uint64_t r = xorshift1024star (rs);
rnumber_4 (dest, kiss); uint32_t high = (uint32_t) (r >> 32);
rnumber_4 (dest, high);
/* Advance to the next element. */ /* Advance to the next element. */
dest += stride0; 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 /* 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 stride0;
index_type dim; index_type dim;
GFC_REAL_8 *dest; GFC_REAL_8 *dest;
GFC_UINTEGER_8 kiss; xorshift1024star_state* rs = get_rand_state();
int n; int n;
dest = x->base_addr; dest = x->base_addr;
@ -449,14 +533,14 @@ arandom_r8 (gfc_array_r8 *x)
stride0 = stride[0]; stride0 = stride[0];
__gthread_mutex_lock (&random_lock); if (unlikely (!rs->init))
init_rand_state (rs, false);
while (dest) while (dest)
{ {
/* random_r8 (dest); */ /* random_r8 (dest); */
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; uint64_t r = xorshift1024star (rs);
kiss += kiss_random_kernel (kiss_seed_2); rnumber_8 (dest, r);
rnumber_8 (dest, kiss);
/* Advance to the next element. */ /* Advance to the next element. */
dest += stride0; dest += stride0;
@ -484,7 +568,6 @@ arandom_r8 (gfc_array_r8 *x)
} }
} }
} }
__gthread_mutex_unlock (&random_lock);
} }
#ifdef HAVE_GFC_REAL_10 #ifdef HAVE_GFC_REAL_10
@ -501,7 +584,7 @@ arandom_r10 (gfc_array_r10 *x)
index_type stride0; index_type stride0;
index_type dim; index_type dim;
GFC_REAL_10 *dest; GFC_REAL_10 *dest;
GFC_UINTEGER_8 kiss; xorshift1024star_state* rs = get_rand_state();
int n; int n;
dest = x->base_addr; dest = x->base_addr;
@ -519,14 +602,14 @@ arandom_r10 (gfc_array_r10 *x)
stride0 = stride[0]; stride0 = stride[0];
__gthread_mutex_lock (&random_lock); if (unlikely (!rs->init))
init_rand_state (rs, false);
while (dest) while (dest)
{ {
/* random_r10 (dest); */ /* random_r10 (dest); */
kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; uint64_t r = xorshift1024star (rs);
kiss += kiss_random_kernel (kiss_seed_2); rnumber_10 (dest, r);
rnumber_10 (dest, kiss);
/* Advance to the next element. */ /* Advance to the next element. */
dest += stride0; dest += stride0;
@ -554,7 +637,6 @@ arandom_r10 (gfc_array_r10 *x)
} }
} }
} }
__gthread_mutex_unlock (&random_lock);
} }
#endif #endif
@ -573,7 +655,7 @@ arandom_r16 (gfc_array_r16 *x)
index_type stride0; index_type stride0;
index_type dim; index_type dim;
GFC_REAL_16 *dest; GFC_REAL_16 *dest;
GFC_UINTEGER_8 kiss1, kiss2; xorshift1024star_state* rs = get_rand_state();
int n; int n;
dest = x->base_addr; dest = x->base_addr;
@ -591,18 +673,15 @@ arandom_r16 (gfc_array_r16 *x)
stride0 = stride[0]; stride0 = stride[0];
__gthread_mutex_lock (&random_lock); if (unlikely (!rs->init))
init_rand_state (rs, false);
while (dest) while (dest)
{ {
/* random_r16 (dest); */ /* random_r16 (dest); */
kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; uint64_t r1 = xorshift1024star (rs);
kiss1 += kiss_random_kernel (kiss_seed_2); uint64_t r2 = xorshift1024star (rs);
rnumber_16 (dest, r1, r2);
kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32;
kiss2 += kiss_random_kernel (kiss_seed_3);
rnumber_16 (dest, kiss1, kiss2);
/* Advance to the next element. */ /* Advance to the next element. */
dest += stride0; dest += stride0;
@ -630,7 +709,6 @@ arandom_r16 (gfc_array_r16 *x)
} }
} }
} }
__gthread_mutex_unlock (&random_lock);
} }
#endif #endif
@ -665,43 +743,17 @@ unscramble_seed (unsigned char *dest, unsigned char *src, int size)
void void
random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
{ {
unsigned char seed[4 * KISS_SIZE]; unsigned char seed[sizeof (master_state)];
#define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_4))
__gthread_mutex_lock (&random_lock);
/* Check that we only have one argument present. */ /* Check that we only have one argument present. */
if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
runtime_error ("RANDOM_SEED should have at most one argument present."); 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];
if (size != NULL) if (size != NULL)
*size = KISS_SIZE; *size = SZ + 1;
if (put != NULL) xorshift1024star_state* rs = get_rand_state();
{
/* 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 the array is too small, abort. */
if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) KISS_SIZE)
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++)
memcpy (seed + i * sizeof(GFC_UINTEGER_4),
&(put->base_addr[(KISS_SIZE - 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. */ /* Return the seed to GET data. */
if (get != NULL) if (get != NULL)
@ -711,20 +763,68 @@ random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
runtime_error ("Array rank of GET is not 1."); runtime_error ("Array rank of GET is not 1.");
/* If the array is too small, abort. */ /* If the array is too small, abort. */
if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) KISS_SIZE) if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1)
runtime_error ("Array size of GET is too small."); runtime_error ("Array size of GET is too small.");
if (!rs->init)
init_rand_state (rs, false);
/* Unscramble the seed. */ /* Unscramble the seed. */
unscramble_seed (seed, (unsigned char *) kiss_seed, 4 * KISS_SIZE); unscramble_seed (seed, (unsigned char *) rs->s, sizeof seed);
/* Then copy it back to the user variable. */ /* Then copy it back to the user variable. */
for (size_t i = 0; i < KISS_SIZE; i++) for (size_t i = 0; i < SZ ; i++)
memcpy (&(get->base_addr[(KISS_SIZE - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]), memcpy (&(get->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]),
seed + i * sizeof(GFC_UINTEGER_4), seed + i * sizeof(GFC_UINTEGER_4),
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)
{
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 (put) != 1)
runtime_error ("Array rank of PUT is not 1.");
/* If the array is too small, abort. */
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 < SZ; i++)
memcpy (seed + i * sizeof(GFC_UINTEGER_4),
&(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 *) master_state, seed, sizeof seed);
njumps = 0;
init_rand_state (rs, true);
rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15;
}
__gthread_mutex_unlock (&random_lock); __gthread_mutex_unlock (&random_lock);
}
#undef SZ
} }
iexport(random_seed_i4); iexport(random_seed_i4);
@ -732,36 +832,15 @@ iexport(random_seed_i4);
void void
random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get) 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. */ /* Check that we only have one argument present. */
if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
runtime_error ("RANDOM_SEED should have at most one argument present."); runtime_error ("RANDOM_SEED should have at most one argument present.");
/* From the standard: "If no argument is present, the processor assigns #define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_8))
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) if (size != NULL)
*size = KISS_SIZE / 2; *size = SZ + 1;
if (put != NULL) xorshift1024star_state* rs = get_rand_state();
{
/* 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 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.");
/* 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)]),
sizeof (GFC_UINTEGER_8));
}
/* Return the seed to GET data. */ /* Return the seed to GET data. */
if (get != NULL) if (get != NULL)
@ -771,24 +850,77 @@ random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
runtime_error ("Array rank of GET is not 1."); runtime_error ("Array rank of GET is not 1.");
/* If the array is too small, abort. */ /* If the array is too small, abort. */
if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) KISS_SIZE / 2) if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1)
runtime_error ("Array size of GET is too small."); runtime_error ("Array size of GET is too small.");
if (!rs->init)
init_rand_state (rs, false);
/* This code now should do correct strides. */ /* This code now should do correct strides. */
for (size_t i = 0; i < KISS_SIZE / 2; i++) for (size_t i = 0; i < SZ; i++)
memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &kiss_seed[2*i], memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &rs->s[i],
sizeof (GFC_UINTEGER_8)); sizeof (GFC_UINTEGER_8));
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)
{
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 (put) != 1)
runtime_error ("Array rank of PUT is not 1.");
/* If the array is too small, abort. */
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 < 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); __gthread_mutex_unlock (&random_lock);
}
} }
iexport(random_seed_i8); iexport(random_seed_i8);
#ifndef __GTHREAD_MUTEX_INIT #if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS
static void __attribute__((constructor)) static void __attribute__((constructor))
init (void) constructor_random (void)
{ {
#ifndef __GTHREAD_MUTEX_INIT
__GTHREAD_MUTEX_INIT_FUNCTION (&random_lock); __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 #endif

2
libgfortran/runtime/main.c
View File

@ -119,8 +119,6 @@ init (void)
set_fpu (); set_fpu ();
init_compile_options (); init_compile_options ();
random_seed_i4 (NULL, NULL, NULL);
} }