123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053 |
- /* Implementation of the RANDOM intrinsics
- Copyright (C) 2002-2022 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).
- Libgfortran is free software; you can redistribute it and/or
- modify it under the terms of the GNU General Public
- License as published by the Free Software Foundation; either
- version 3 of the License, or (at your option) any later version.
- Ligbfortran is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
- Under Section 7 of GPL version 3, you are granted additional
- permissions described in the GCC Runtime Library Exception, version
- 3.1, as published by the Free Software Foundation.
- You should have received a copy of the GNU General Public License and
- a copy of the GCC Runtime Library Exception along with this program;
- 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);
- extern void random_r8 (GFC_REAL_8 *);
- iexport_proto(random_r8);
- extern void arandom_r4 (gfc_array_r4 *);
- export_proto(arandom_r4);
- extern void arandom_r8 (gfc_array_r8 *);
- export_proto(arandom_r8);
- #ifdef HAVE_GFC_REAL_10
- extern void random_r10 (GFC_REAL_10 *);
- iexport_proto(random_r10);
- extern void arandom_r10 (gfc_array_r10 *);
- export_proto(arandom_r10);
- #endif
- #ifdef HAVE_GFC_REAL_16
- extern void random_r16 (GFC_REAL_16 *);
- iexport_proto(random_r16);
- extern void arandom_r16 (gfc_array_r16 *);
- export_proto(arandom_r16);
- #endif
- #ifdef HAVE_GFC_REAL_17
- extern void random_r17 (GFC_REAL_17 *);
- iexport_proto(random_r17);
- extern void arandom_r17 (gfc_array_r17 *);
- export_proto(arandom_r17);
- #endif
- #ifdef __GTHREAD_MUTEX_INIT
- static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
- #else
- static __gthread_mutex_t random_lock;
- #endif
- /* Helper routines to map a GFC_UINTEGER_* to the corresponding
- GFC_REAL_* types in the range of [0,1). If GFC_REAL_*_RADIX are 2
- or 16, respectively, we mask off the bits that don't fit into the
- correct GFC_REAL_*, convert to the real type, then multiply by the
- correct offset. */
- static void
- rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
- {
- GFC_UINTEGER_4 mask;
- #if GFC_REAL_4_RADIX == 2
- mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS);
- #elif GFC_REAL_4_RADIX == 16
- mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4);
- #else
- #error "GFC_REAL_4_RADIX has unknown value"
- #endif
- v = v & mask;
- *f = (GFC_REAL_4) v * GFC_REAL_4_LITERAL(0x1.p-32);
- }
- static void
- rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
- {
- GFC_UINTEGER_8 mask;
- #if GFC_REAL_8_RADIX == 2
- mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS);
- #elif GFC_REAL_8_RADIX == 16
- mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4);
- #else
- #error "GFC_REAL_8_RADIX has unknown value"
- #endif
- v = v & mask;
- *f = (GFC_REAL_8) v * GFC_REAL_8_LITERAL(0x1.p-64);
- }
- #ifdef HAVE_GFC_REAL_10
- static void
- rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
- {
- GFC_UINTEGER_8 mask;
- #if GFC_REAL_10_RADIX == 2
- mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS);
- #elif GFC_REAL_10_RADIX == 16
- mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4);
- #else
- #error "GFC_REAL_10_RADIX has unknown value"
- #endif
- v = v & mask;
- *f = (GFC_REAL_10) v * GFC_REAL_10_LITERAL(0x1.p-64);
- }
- #endif
- #ifdef HAVE_GFC_REAL_16
- /* For REAL(KIND=16), we only need to mask off the lower bits. */
- static void
- rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
- {
- GFC_UINTEGER_8 mask;
- #if GFC_REAL_16_RADIX == 2
- mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS);
- #elif GFC_REAL_16_RADIX == 16
- mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4);
- #else
- #error "GFC_REAL_16_RADIX has unknown value"
- #endif
- v2 = v2 & mask;
- *f = (GFC_REAL_16) v1 * GFC_REAL_16_LITERAL(0x1.p-64)
- + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(0x1.p-128);
- }
- #endif
- #ifdef HAVE_GFC_REAL_17
- /* For REAL(KIND=16), we only need to mask off the lower bits. */
- static void
- rnumber_17 (GFC_REAL_17 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
- {
- GFC_UINTEGER_8 mask;
- #if GFC_REAL_17_RADIX == 2
- mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_17_DIGITS);
- #elif GFC_REAL_17_RADIX == 16
- mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_17_DIGITS) * 4);
- #else
- #error "GFC_REAL_17_RADIX has unknown value"
- #endif
- v2 = v2 & mask;
- *f = (GFC_REAL_17) v1 * GFC_REAL_17_LITERAL(0x1.p-64)
- + (GFC_REAL_17) v2 * GFC_REAL_17_LITERAL(0x1.p-128);
- }
- #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 __gthread_key_t rand_state_key;
- 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));
- }
- static uint64_t
- prng_next (prng_state* rs)
- {
- const uint64_t result = rotl(rs->s[1] * 5, 7) * 9;
- const uint64_t t = rs->s[1] << 17;
- rs->s[2] ^= rs->s[0];
- rs->s[3] ^= rs->s[1];
- rs->s[1] ^= rs->s[2];
- rs->s[0] ^= rs->s[3];
- 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)
- {
- 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);
- /* This function produces a REAL(8) value from the uniform distribution
- with range [0,1). */
- void
- random_r8 (GFC_REAL_8 *x)
- {
- GFC_UINTEGER_8 r;
- prng_state* rs = get_rand_state();
- if (unlikely (!rs->init))
- init_rand_state (rs, false);
- r = prng_next (rs);
- rnumber_8 (x, r);
- }
- iexport(random_r8);
- #ifdef HAVE_GFC_REAL_10
- /* This function produces a REAL(10) value from the uniform distribution
- with range [0,1). */
- void
- random_r10 (GFC_REAL_10 *x)
- {
- GFC_UINTEGER_8 r;
- prng_state* rs = get_rand_state();
- if (unlikely (!rs->init))
- init_rand_state (rs, false);
- r = prng_next (rs);
- rnumber_10 (x, r);
- }
- iexport(random_r10);
- #endif
- /* This function produces a REAL(16) value from the uniform distribution
- with range [0,1). */
- #ifdef HAVE_GFC_REAL_16
- void
- random_r16 (GFC_REAL_16 *x)
- {
- 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 produces a REAL(16) value from the uniform distribution
- with range [0,1). */
- #ifdef HAVE_GFC_REAL_17
- void
- random_r17 (GFC_REAL_17 *x)
- {
- 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_17 (x, r1, r2);
- }
- iexport(random_r17);
- #endif
- /* This function fills a REAL(4) array with values from the uniform
- distribution with range [0,1). */
- void
- arandom_r4 (gfc_array_r4 *x)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type stride[GFC_MAX_DIMENSIONS];
- index_type stride0;
- index_type dim;
- GFC_REAL_4 *dest;
- prng_state* rs = get_rand_state();
- dest = x->base_addr;
- dim = GFC_DESCRIPTOR_RANK (x);
- for (index_type n = 0; n < dim; n++)
- {
- count[n] = 0;
- stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
- extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
- if (extent[n] <= 0)
- return;
- }
- stride0 = stride[0];
- if (unlikely (!rs->init))
- init_rand_state (rs, false);
- while (dest)
- {
- /* random_r4 (dest); */
- 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. */
- index_type n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- dest -= stride[n] * extent[n];
- n++;
- if (n == dim)
- {
- dest = NULL;
- break;
- }
- else
- {
- count[n]++;
- dest += stride[n];
- }
- }
- }
- }
- /* This function fills a REAL(8) array with values from the uniform
- distribution with range [0,1). */
- void
- arandom_r8 (gfc_array_r8 *x)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type stride[GFC_MAX_DIMENSIONS];
- index_type stride0;
- index_type dim;
- GFC_REAL_8 *dest;
- prng_state* rs = get_rand_state();
- dest = x->base_addr;
- dim = GFC_DESCRIPTOR_RANK (x);
- for (index_type n = 0; n < dim; n++)
- {
- count[n] = 0;
- stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
- extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
- if (extent[n] <= 0)
- return;
- }
- stride0 = stride[0];
- if (unlikely (!rs->init))
- init_rand_state (rs, false);
- while (dest)
- {
- /* random_r8 (dest); */
- 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. */
- index_type n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- dest -= stride[n] * extent[n];
- n++;
- if (n == dim)
- {
- dest = NULL;
- break;
- }
- else
- {
- count[n]++;
- dest += stride[n];
- }
- }
- }
- }
- #ifdef HAVE_GFC_REAL_10
- /* This function fills a REAL(10) array with values from the uniform
- distribution with range [0,1). */
- void
- arandom_r10 (gfc_array_r10 *x)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type stride[GFC_MAX_DIMENSIONS];
- index_type stride0;
- index_type dim;
- GFC_REAL_10 *dest;
- prng_state* rs = get_rand_state();
- dest = x->base_addr;
- dim = GFC_DESCRIPTOR_RANK (x);
- for (index_type n = 0; n < dim; n++)
- {
- count[n] = 0;
- stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
- extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
- if (extent[n] <= 0)
- return;
- }
- stride0 = stride[0];
- if (unlikely (!rs->init))
- init_rand_state (rs, false);
- while (dest)
- {
- /* random_r10 (dest); */
- 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. */
- index_type n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- dest -= stride[n] * extent[n];
- n++;
- if (n == dim)
- {
- dest = NULL;
- break;
- }
- else
- {
- count[n]++;
- dest += stride[n];
- }
- }
- }
- }
- #endif
- #ifdef HAVE_GFC_REAL_16
- /* This function fills a REAL(16) array with values from the uniform
- distribution with range [0,1). */
- void
- arandom_r16 (gfc_array_r16 *x)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type stride[GFC_MAX_DIMENSIONS];
- index_type stride0;
- index_type dim;
- GFC_REAL_16 *dest;
- prng_state* rs = get_rand_state();
- dest = x->base_addr;
- dim = GFC_DESCRIPTOR_RANK (x);
- for (index_type n = 0; n < dim; n++)
- {
- count[n] = 0;
- stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
- extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
- if (extent[n] <= 0)
- return;
- }
- stride0 = stride[0];
- if (unlikely (!rs->init))
- init_rand_state (rs, false);
- while (dest)
- {
- /* random_r16 (dest); */
- 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. */
- index_type n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- dest -= stride[n] * extent[n];
- n++;
- if (n == dim)
- {
- dest = NULL;
- break;
- }
- else
- {
- count[n]++;
- dest += stride[n];
- }
- }
- }
- }
- #endif
- #ifdef HAVE_GFC_REAL_17
- /* This function fills a REAL(16) array with values from the uniform
- distribution with range [0,1). */
- void
- arandom_r17 (gfc_array_r17 *x)
- {
- index_type count[GFC_MAX_DIMENSIONS];
- index_type extent[GFC_MAX_DIMENSIONS];
- index_type stride[GFC_MAX_DIMENSIONS];
- index_type stride0;
- index_type dim;
- GFC_REAL_17 *dest;
- prng_state* rs = get_rand_state();
- dest = x->base_addr;
- dim = GFC_DESCRIPTOR_RANK (x);
- for (index_type n = 0; n < dim; n++)
- {
- count[n] = 0;
- stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
- extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
- if (extent[n] <= 0)
- return;
- }
- stride0 = stride[0];
- if (unlikely (!rs->init))
- init_rand_state (rs, false);
- while (dest)
- {
- /* random_r17 (dest); */
- uint64_t r1 = prng_next (rs);
- uint64_t r2 = prng_next (rs);
- rnumber_17 (dest, r1, r2);
- /* Advance to the next element. */
- dest += stride0;
- count[0]++;
- /* Advance to the next source element. */
- index_type n = 0;
- while (count[n] == extent[n])
- {
- /* When we get to the end of a dimension, reset it and increment
- the next dimension. */
- count[n] = 0;
- /* We could precalculate these products, but this is a less
- frequently used path so probably not worth it. */
- dest -= stride[n] * extent[n];
- n++;
- if (n == dim)
- {
- dest = NULL;
- break;
- }
- else
- {
- count[n]++;
- dest += stride[n];
- }
- }
- }
- }
- #endif
- /* Number of elements in master_state array. */
- #define SZU64 (sizeof (master_state.s) / sizeof (uint64_t))
- /* Equivalent number of elements in an array of GFC_INTEGER_{4,8}. */
- #define SZ_IN_INT_4 (SZU64 * (sizeof (uint64_t) / sizeof (GFC_INTEGER_4)))
- #define SZ_IN_INT_8 (SZU64 * (sizeof (uint64_t) / sizeof (GFC_INTEGER_8)))
- /* 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
- scramble_seed (uint64_t *dest, const uint64_t *src)
- {
- 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)
- {
- 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.");
- if (size != NULL)
- *size = SZ_IN_INT_4;
- 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_IN_INT_4)
- 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_IN_INT_4 ; i++)
- memcpy (&(get->base_addr[(SZ_IN_INT_4 - 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)
- {
- master_state.init = false;
- 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_IN_INT_4)
- runtime_error ("Array size of PUT is too small.");
- /* We copy the seed given by the user. */
- for (size_t i = 0; i < SZ_IN_INT_4; i++)
- memcpy ((unsigned char*) seed + i * sizeof(GFC_UINTEGER_4),
- &(put->base_addr[(SZ_IN_INT_4 - 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 (master_state.s, seed);
- master_state.init = true;
- init_rand_state (rs, true);
- }
- __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)
- {
- 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.");
- if (size != NULL)
- *size = SZ_IN_INT_8;
- 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_IN_INT_8)
- 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);
- /* This code now should do correct strides. */
- for (size_t i = 0; i < SZ_IN_INT_8; i++)
- memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &seed[i],
- sizeof (GFC_UINTEGER_8));
- }
- 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)
- {
- master_state.init = false;
- 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_IN_INT_8)
- runtime_error ("Array size of PUT is too small.");
- /* This code now should do correct strides. */
- for (size_t i = 0; i < SZ_IN_INT_8; i++)
- memcpy (&seed[i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]),
- sizeof (GFC_UINTEGER_8));
- scramble_seed (master_state.s, seed);
- master_state.init = true;
- init_rand_state (rs, true);
- }
- __gthread_mutex_unlock (&random_lock);
- }
- }
- iexport(random_seed_i8);
- #if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS
- static void __attribute__((constructor))
- 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
|