/* Implementation of the RANDOM intrinsics
- Copyright (C) 2002-2013 Free Software Foundation, Inc.
- Contributed by Lars Segerlund <seger@linuxmail.org>
- and Steve Kargl.
+ Copyright (C) 2002-2020 Free Software Foundation, Inc.
+ Contributed by Lars Segerlund <seger@linuxmail.org>,
+ Steve Kargl and Janne Blomqvist.
This file is part of the GNU Fortran runtime library (libgfortran).
see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
<http://www.gnu.org/licenses/>. */
+/* For rand_s. */
+#define _CRT_RAND_S
+
#include "libgfortran.h"
#include <gthr.h>
#include <string.h>
+#ifdef HAVE_UNISTD_H
+#include <unistd.h>
+#endif
+#include <sys/stat.h>
+#include <fcntl.h>
+#include "time_1.h"
+#ifdef HAVE_SYS_RANDOM_H
+#include <sys/random.h>
+#endif
+
+#ifdef __MINGW32__
+#define HAVE_GETPID 1
+#include <process.h>
+#include <_mingw.h> /* For __MINGW64_VERSION_MAJOR */
+#endif
+
extern void random_r4 (GFC_REAL_4 *);
iexport_proto(random_r4);
+ (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 <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
- 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" <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
+/*
+
+ We use the xoshiro256** generator, a fast high-quality generator
+ that:
+
+ - passes TestU1 without any failures
+
+ - provides a "jump" function making it easy to provide many
+ independent parallel streams.
+
+ - Long period of 2**256 - 1
+
+ A description can be found at
+
+ http://prng.di.unimi.it/
+
+ or
+
+ https://arxiv.org/abs/1805.01407
+
+ The paper includes public domain source code which is the basis for
+ the implementation below.
+
+*/
+typedef struct
+{
+ bool init;
+ uint64_t s[4];
+}
+prng_state;
+
+
+/* master_state is the only variable protected by random_lock. */
+static prng_state master_state = { .init = false, .s = {
+ 0xad63fa1ed3b55f36ULL, 0xd94473e78978b497ULL, 0xbc60592a98172477ULL,
+ 0xa3de7c6e81265301ULL }
};
-static const GFC_INTEGER_4 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;
+static __gthread_key_t rand_state_key;
-#ifdef HAVE_GFC_REAL_16
-static GFC_UINTEGER_4 * const kiss_seed_3 = kiss_seed + 8;
-#endif
+static prng_state*
+get_rand_state (void)
+{
+ /* For single threaded apps. */
+ static prng_state rand_state;
+
+ if (__gthread_active_p ())
+ {
+ void* p = __gthread_getspecific (rand_state_key);
+ if (!p)
+ {
+ p = xcalloc (1, sizeof (prng_state));
+ __gthread_setspecific (rand_state_key, p);
+ }
+ return p;
+ }
+ else
+ return &rand_state;
+}
+
+static inline uint64_t
+rotl (const uint64_t x, int k)
+{
+ return (x << k) | (x >> (64 - k));
+}
-/* 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)
+static uint64_t
+prng_next (prng_state* rs)
{
- GFC_UINTEGER_4 kiss;
+ const uint64_t result = rotl(rs->s[1] * 5, 7) * 9;
+
+ const uint64_t t = rs->s[1] << 17;
- 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];
+ rs->s[2] ^= rs->s[0];
+ rs->s[3] ^= rs->s[1];
+ rs->s[1] ^= rs->s[2];
+ rs->s[0] ^= rs->s[3];
- return kiss;
+ rs->s[2] ^= t;
+
+ rs->s[3] = rotl(rs->s[3], 45);
+
+ return result;
+}
+
+
+/* This is the jump function for the generator. It is equivalent to
+ 2^128 calls to prng_next(); it can be used to generate 2^128
+ non-overlapping subsequences for parallel computations. */
+
+static void
+jump (prng_state* rs)
+{
+ static const uint64_t JUMP[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c };
+
+ uint64_t s0 = 0;
+ uint64_t s1 = 0;
+ uint64_t s2 = 0;
+ uint64_t s3 = 0;
+ for(size_t i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
+ for(int b = 0; b < 64; b++) {
+ if (JUMP[i] & UINT64_C(1) << b) {
+ s0 ^= rs->s[0];
+ s1 ^= rs->s[1];
+ s2 ^= rs->s[2];
+ s3 ^= rs->s[3];
+ }
+ prng_next (rs);
+ }
+
+ rs->s[0] = s0;
+ rs->s[1] = s1;
+ rs->s[2] = s2;
+ rs->s[3] = s3;
}
+
+/* Splitmix64 recommended by xoshiro author for initializing. After
+ getting one uint64_t value from the OS, this is used to fill in the
+ rest of the xoshiro state. */
+
+static uint64_t
+splitmix64 (uint64_t x)
+{
+ uint64_t z = (x += 0x9e3779b97f4a7c15);
+ z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
+ z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
+ return z ^ (z >> 31);
+}
+
+
+/* Get some bytes from the operating system in order to seed
+ the PRNG. */
+
+static int
+getosrandom (void *buf, size_t buflen)
+{
+ /* rand_s is available in MinGW-w64 but not plain MinGW. */
+#if defined(__MINGW64_VERSION_MAJOR)
+ unsigned int* b = buf;
+ for (size_t i = 0; i < buflen / sizeof (unsigned int); i++)
+ rand_s (&b[i]);
+ return buflen;
+#else
+#ifdef HAVE_GETENTROPY
+ if (getentropy (buf, buflen) == 0)
+ return buflen;
+#endif
+ 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;
+ }
+ uint64_t seed = 0x047f7684e9fc949dULL;
+ 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
+ size_t size = buflen < sizeof (uint64_t) ? buflen : sizeof (uint64_t);
+ memcpy (buf, &seed, size);
+ return size;
+#endif /* __MINGW64_VERSION_MAJOR */
+}
+
+
+/* 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 (prng_state* rs, const bool locked)
+{
+ if (!locked)
+ __gthread_mutex_lock (&random_lock);
+ if (!master_state.init)
+ {
+ uint64_t os_seed;
+ getosrandom (&os_seed, sizeof (os_seed));
+ for (uint64_t i = 0; i < sizeof (master_state.s) / sizeof (uint64_t); i++)
+ {
+ os_seed = splitmix64 (os_seed);
+ master_state.s[i] = os_seed;
+ }
+ master_state.init = true;
+ }
+ memcpy (&rs->s, master_state.s, sizeof (master_state.s));
+ jump (&master_state);
+ if (!locked)
+ __gthread_mutex_unlock (&random_lock);
+ rs->init = true;
+}
+
+
/* 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);
+ prng_state* rs = get_rand_state();
+
+ if (unlikely (!rs->init))
+ init_rand_state (rs, false);
+ uint64_t r = prng_next (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);
void
random_r8 (GFC_REAL_8 *x)
{
- GFC_UINTEGER_8 kiss;
+ GFC_UINTEGER_8 r;
+ prng_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 = prng_next (rs);
+ rnumber_8 (x, r);
}
iexport(random_r8);
void
random_r10 (GFC_REAL_10 *x)
{
- GFC_UINTEGER_8 kiss;
+ GFC_UINTEGER_8 r;
+ prng_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 = prng_next (rs);
+ rnumber_10 (x, r);
}
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;
+ prng_state* rs = get_rand_state();
+
+ if (unlikely (!rs->init))
+ init_rand_state (rs, false);
+ r1 = prng_next (rs);
+ r2 = prng_next (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). */
index_type stride0;
index_type dim;
GFC_REAL_4 *dest;
- GFC_UINTEGER_4 kiss;
- int n;
+ prng_state* rs = get_rand_state();
dest = x->base_addr;
dim = GFC_DESCRIPTOR_RANK (x);
- for (n = 0; n < dim; n++)
+ for (index_type n = 0; n < dim; n++)
{
count[n] = 0;
stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
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 = prng_next (rs);
+ uint32_t high = (uint32_t) (r >> 32);
+ rnumber_4 (dest, high);
/* Advance to the next element. */
dest += stride0;
count[0]++;
/* Advance to the next source element. */
- n = 0;
+ index_type n = 0;
while (count[n] == extent[n])
{
/* When we get to the end of a dimension, reset it and increment
}
}
}
- __gthread_mutex_unlock (&random_lock);
}
/* This function fills a REAL(8) array with values from the uniform
index_type stride0;
index_type dim;
GFC_REAL_8 *dest;
- GFC_UINTEGER_8 kiss;
- int n;
+ prng_state* rs = get_rand_state();
dest = x->base_addr;
dim = GFC_DESCRIPTOR_RANK (x);
- for (n = 0; n < dim; n++)
+ for (index_type n = 0; n < dim; n++)
{
count[n] = 0;
stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
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 = prng_next (rs);
+ rnumber_8 (dest, r);
/* Advance to the next element. */
dest += stride0;
count[0]++;
/* Advance to the next source element. */
- n = 0;
+ index_type n = 0;
while (count[n] == extent[n])
{
/* When we get to the end of a dimension, reset it and increment
}
}
}
- __gthread_mutex_unlock (&random_lock);
}
#ifdef HAVE_GFC_REAL_10
index_type stride0;
index_type dim;
GFC_REAL_10 *dest;
- GFC_UINTEGER_8 kiss;
- int n;
+ prng_state* rs = get_rand_state();
dest = x->base_addr;
dim = GFC_DESCRIPTOR_RANK (x);
- for (n = 0; n < dim; n++)
+ for (index_type n = 0; n < dim; n++)
{
count[n] = 0;
stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
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 = prng_next (rs);
+ rnumber_10 (dest, r);
/* Advance to the next element. */
dest += stride0;
count[0]++;
/* Advance to the next source element. */
- n = 0;
+ index_type n = 0;
while (count[n] == extent[n])
{
/* When we get to the end of a dimension, reset it and increment
}
}
}
- __gthread_mutex_unlock (&random_lock);
}
#endif
index_type stride0;
index_type dim;
GFC_REAL_16 *dest;
- GFC_UINTEGER_8 kiss1, kiss2;
- int n;
+ prng_state* rs = get_rand_state();
dest = x->base_addr;
dim = GFC_DESCRIPTOR_RANK (x);
- for (n = 0; n < dim; n++)
+ for (index_type n = 0; n < dim; n++)
{
count[n] = 0;
stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
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 = prng_next (rs);
+ uint64_t r2 = prng_next (rs);
+ rnumber_16 (dest, r1, r2);
/* Advance to the next element. */
dest += stride0;
count[0]++;
/* Advance to the next source element. */
- n = 0;
+ index_type n = 0;
while (count[n] == extent[n])
{
/* When we get to the end of a dimension, reset it and increment
}
}
}
- __gthread_mutex_unlock (&random_lock);
}
#endif
+/* Number of elements in master_state array. */
+#define SZU64 (sizeof (master_state.s) / sizeof (uint64_t))
-static void
-scramble_seed (unsigned char *dest, unsigned char *src, int size)
-{
- int i;
- for (i = 0; i < size; i++)
- dest[(i % 2) * (size / 2) + i / 2] = src[i];
-}
+/* Keys for scrambling the seed in order to avoid poor seeds. */
+static const uint64_t xor_keys[] = {
+ 0xbd0c5b6e50c2df49ULL, 0xd46061cd46e1df38ULL, 0xbb4f4d4ed6103544ULL,
+ 0x114a583d0756ad39ULL
+};
+
+
+/* Since a XOR cipher is symmetric, we need only one routine, and we
+ can use it both for encryption and decryption. */
static void
-unscramble_seed (unsigned char *dest, unsigned char *src, int size)
+scramble_seed (uint64_t *dest, const uint64_t *src)
{
- int i;
-
- for (i = 0; i < size; i++)
- dest[i] = src[(i % 2) * (size / 2) + i / 2];
+ for (size_t i = 0; i < SZU64; i++)
+ dest[i] = src[i] ^ xor_keys[i];
}
-
/* random_seed is used to seed the PRNG with either a default
set of seeds or user specified set of seeds. random_seed
must be called with no argument or exactly one argument. */
void
random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
{
- int i;
- unsigned char seed[4*kiss_size];
-
- __gthread_mutex_lock (&random_lock);
+ uint64_t seed[SZU64];
+#define SZ (sizeof (master_state.s) / 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;
+
+ prng_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)
+ runtime_error ("Array size of GET is too small.");
+
+ if (!rs->init)
+ init_rand_state (rs, false);
+
+ /* Unscramble the seed. */
+ scramble_seed (seed, rs->s);
+
+ /* 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)]),
+ (unsigned char*) seed + i * sizeof(GFC_UINTEGER_4),
+ sizeof(GFC_UINTEGER_4));
+ }
+
+ 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 (i = 0; i < kiss_size; i++)
- kiss_seed[i] = kiss_default_seed[i];
-
- if (size != NULL)
- *size = kiss_size;
+ {
+ master_state.init = false;
+ init_rand_state (rs, true);
+ }
if (put != NULL)
{
runtime_error ("Array rank of PUT is not 1.");
/* If the array is too small, abort. */
- if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size)
+ if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ)
runtime_error ("Array size of PUT is too small.");
/* We copy the seed given by the user. */
- for (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)]),
+ for (size_t i = 0; i < SZ; i++)
+ memcpy ((unsigned char*) 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 *) kiss_seed, seed, 4*kiss_size);
+ scramble_seed (master_state.s, seed);
+ master_state.init = true;
+ init_rand_state (rs, true);
}
+ __gthread_mutex_unlock (&random_lock);
+ }
+#undef SZ
+}
+iexport(random_seed_i4);
+
+
+void
+random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
+{
+ uint64_t seed[SZU64];
+
+ /* 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.");
+
+#define SZ (sizeof (master_state.s) / sizeof (GFC_INTEGER_8))
+ if (size != NULL)
+ *size = SZ;
+
+ prng_state* rs = get_rand_state();
+
/* Return the seed to GET data. */
if (get != NULL)
{
runtime_error ("Array rank of GET is not 1.");
/* If the array is too small, abort. */
- if (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size)
+ if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ)
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 *) kiss_seed, 4*kiss_size);
+ scramble_seed (seed, rs->s);
- /* Then copy it back to the user variable. */
- for (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));
+ /* This code now should do correct strides. */
+ for (size_t i = 0; i < SZ; i++)
+ memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &seed[i],
+ sizeof (GFC_UINTEGER_8));
}
- __gthread_mutex_unlock (&random_lock);
-}
-iexport(random_seed_i4);
-
-
-void
-random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
-{
- int i;
-
+ else
+ {
__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 (i = 0; i < kiss_size; i++)
- kiss_seed[i] = kiss_default_seed[i];
-
- if (size != NULL)
- *size = kiss_size / 2;
+ {
+ master_state.init = false;
+ init_rand_state (rs, true);
+ }
if (put != NULL)
{
runtime_error ("Array rank of PUT is not 1.");
/* If the array is too small, abort. */
- if (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size / 2)
+ if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ)
runtime_error ("Array size of PUT is too small.");
/* This code now should do correct strides. */
- for (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 (&seed[i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]),
sizeof (GFC_UINTEGER_8));
- }
-
- /* 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) < kiss_size / 2)
- runtime_error ("Array size of GET is too small.");
+ scramble_seed (master_state.s, seed);
+ master_state.init = true;
+ init_rand_state (rs, true);
+ }
- /* This code now should do correct strides. */
- for (i = 0; i < kiss_size / 2; i++)
- memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &kiss_seed[2*i],
- sizeof (GFC_UINTEGER_8));
- }
__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
projects / thirdparty / gcc.git / blobdiff