libgfortran/intrinsics/random.c

   1 /* Implementation of the RANDOM intrinsics
   2    Copyright (C) 2002-2019 Free Software Foundation, Inc.
   3    Contributed by Lars Segerlund <seger@linuxmail.org>,
   4    Steve Kargl and Janne Blomqvist.
   5
   6 This file is part of the GNU Fortran runtime library (libgfortran).
   7
   8 Libgfortran is free software; you can redistribute it and/or
   9 modify it under the terms of the GNU General Public
  10 License as published by the Free Software Foundation; either
  11 version 3 of the License, or (at your option) any later version.
  12
  13 Ligbfortran is distributed in the hope that it will be useful,
  14 but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16 GNU General Public License for more details.
  17
  18 Under Section 7 of GPL version 3, you are granted additional
  19 permissions described in the GCC Runtime Library Exception, version
  20 3.1, as published by the Free Software Foundation.
  21
  22 You should have received a copy of the GNU General Public License and
  23 a copy of the GCC Runtime Library Exception along with this program;
  24 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  25 <http://www.gnu.org/licenses/>.  */
  26
  27 /* For rand_s.  */
  28 #define _CRT_RAND_S
  29
  30 #include "libgfortran.h"
  31 #include <gthr.h>
  32 #include <string.h>
  33
  34 #ifdef HAVE_UNISTD_H
  35 #include <unistd.h>
  36 #endif
  37 #include <sys/stat.h>
  38 #include <fcntl.h>
  39 #include "time_1.h"
  40 #ifdef HAVE_SYS_RANDOM_H
  41 #include <sys/random.h>
  42 #endif
  43
  44 #ifdef __MINGW32__
  45 #define HAVE_GETPID 1
  46 #include <process.h>
  47 #include <_mingw.h> /* For __MINGW64_VERSION_MAJOR  */
  48 #endif
  49
  50 extern void random_r4 (GFC_REAL_4 *);
  51 iexport_proto(random_r4);
  52
  53 extern void random_r8 (GFC_REAL_8 *);
  54 iexport_proto(random_r8);
  55
  56 extern void arandom_r4 (gfc_array_r4 *);
  57 export_proto(arandom_r4);
  58
  59 extern void arandom_r8 (gfc_array_r8 *);
  60 export_proto(arandom_r8);
  61
  62 #ifdef HAVE_GFC_REAL_10
  63
  64 extern void random_r10 (GFC_REAL_10 *);
  65 iexport_proto(random_r10);
  66
  67 extern void arandom_r10 (gfc_array_r10 *);
  68 export_proto(arandom_r10);
  69
  70 #endif
  71
  72 #ifdef HAVE_GFC_REAL_16
  73
  74 extern void random_r16 (GFC_REAL_16 *);
  75 iexport_proto(random_r16);
  76
  77 extern void arandom_r16 (gfc_array_r16 *);
  78 export_proto(arandom_r16);
  79
  80 #endif
  81
  82 #ifdef __GTHREAD_MUTEX_INIT
  83 static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT;
  84 #else
  85 static __gthread_mutex_t random_lock;
  86 #endif
  87
  88 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
  89    GFC_REAL_* types in the range of [0,1).  If GFC_REAL_*_RADIX are 2
  90    or 16, respectively, we mask off the bits that don't fit into the
  91    correct GFC_REAL_*, convert to the real type, then multiply by the
  92    correct offset.  */
  93
  94
  95 static void
  96 rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v)
  97 {
  98   GFC_UINTEGER_4 mask;
  99 #if GFC_REAL_4_RADIX == 2
 100   mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS);
 101 #elif GFC_REAL_4_RADIX == 16
 102   mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4);
 103 #else
 104 #error "GFC_REAL_4_RADIX has unknown value"
 105 #endif
 106   v = v & mask;
 107   *f = (GFC_REAL_4) v * GFC_REAL_4_LITERAL(0x1.p-32);
 108 }
 109
 110 static void
 111 rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v)
 112 {
 113   GFC_UINTEGER_8 mask;
 114 #if GFC_REAL_8_RADIX == 2
 115   mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS);
 116 #elif GFC_REAL_8_RADIX == 16
 117   mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4);
 118 #else
 119 #error "GFC_REAL_8_RADIX has unknown value"
 120 #endif
 121   v = v & mask;
 122   *f = (GFC_REAL_8) v * GFC_REAL_8_LITERAL(0x1.p-64);
 123 }
 124
 125 #ifdef HAVE_GFC_REAL_10
 126
 127 static void
 128 rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v)
 129 {
 130   GFC_UINTEGER_8 mask;
 131 #if GFC_REAL_10_RADIX == 2
 132   mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS);
 133 #elif GFC_REAL_10_RADIX == 16
 134   mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4);
 135 #else
 136 #error "GFC_REAL_10_RADIX has unknown value"
 137 #endif
 138   v = v & mask;
 139   *f = (GFC_REAL_10) v * GFC_REAL_10_LITERAL(0x1.p-64);
 140 }
 141 #endif
 142
 143 #ifdef HAVE_GFC_REAL_16
 144
 145 /* For REAL(KIND=16), we only need to mask off the lower bits.  */
 146
 147 static void
 148 rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2)
 149 {
 150   GFC_UINTEGER_8 mask;
 151 #if GFC_REAL_16_RADIX == 2
 152   mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS);
 153 #elif GFC_REAL_16_RADIX == 16
 154   mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4);
 155 #else
 156 #error "GFC_REAL_16_RADIX has unknown value"
 157 #endif
 158   v2 = v2 & mask;
 159   *f = (GFC_REAL_16) v1 * GFC_REAL_16_LITERAL(0x1.p-64)
 160     + (GFC_REAL_16) v2 * GFC_REAL_16_LITERAL(0x1.p-128);
 161 }
 162 #endif
 163
 164
 165 /*
 166
 167    We use the xorshift1024* generator, a fast high-quality generator
 168    that:
 169
 170    - passes TestU1 without any failures
 171
 172    - provides a "jump" function making it easy to provide many
 173      independent parallel streams.
 174
 175    - Long period of 2**1024 - 1
 176
 177    A description can be found at
 178
 179    http://vigna.di.unimi.it/ftp/papers/xorshift.pdf
 180
 181    or
 182
 183    http://arxiv.org/abs/1402.6246
 184
 185    The paper includes public domain source code which is the basis for
 186    the implementation below.
 187
 188 */
 189 typedef struct
 190 {
 191   bool init;
 192   int p;
 193   uint64_t s[16];
 194 }
 195 xorshift1024star_state;
 196
 197
 198 /* master_init, njumps, and master_state are the only variables
 199    protected by random_lock.  */
 200 static bool master_init;
 201 static unsigned njumps; /* How many times we have jumped.  */
 202 static uint64_t master_state[] = {
 203   0xad63fa1ed3b55f36ULL, 0xd94473e78978b497ULL, 0xbc60592a98172477ULL,
 204   0xa3de7c6e81265301ULL, 0x586640c5e785af27ULL, 0x7a2a3f63b67ce5eaULL,
 205   0x9fde969f922d9b82ULL, 0xe6fe34379b3f3822ULL, 0x6c277eac3e99b6c2ULL,
 206   0x9197290ab0d3f069ULL, 0xdb227302f6c25576ULL, 0xee0209aee527fae9ULL,
 207   0x675666a793cd05b9ULL, 0xd048c99fbc70c20fULL, 0x775f8c3dba385ef5ULL,
 208   0x625288bc262faf33ULL
 209 };
 210
 211
 212 static __gthread_key_t rand_state_key;
 213
 214 static xorshift1024star_state*
 215 get_rand_state (void)
 216 {
 217   /* For single threaded apps.  */
 218   static xorshift1024star_state rand_state;
 219
 220   if (__gthread_active_p ())
 221     {
 222       void* p = __gthread_getspecific (rand_state_key);
 223       if (!p)
 224         {
 225           p = xcalloc (1, sizeof (xorshift1024star_state));
 226           __gthread_setspecific (rand_state_key, p);
 227         }
 228       return p;
 229     }
 230   else
 231     return &rand_state;
 232 }
 233
 234
 235 static uint64_t
 236 xorshift1024star (xorshift1024star_state* rs)
 237 {
 238   int p = rs->p;
 239   const uint64_t s0 = rs->s[p];
 240   uint64_t s1 = rs->s[p = (p + 1) & 15];
 241   s1 ^= s1 << 31;
 242   rs->s[p] = s1 ^ s0 ^ (s1 >> 11) ^ (s0 >> 30);
 243   rs->p = p;
 244   return rs->s[p] * UINT64_C(1181783497276652981);
 245 }
 246
 247
 248 /* This is the jump function for the generator. It is equivalent to
 249    2^512 calls to xorshift1024star(); it can be used to generate 2^512
 250    non-overlapping subsequences for parallel computations. */
 251
 252 static void
 253 jump (xorshift1024star_state* rs)
 254 {
 255   static const uint64_t JUMP[] = {
 256     0x84242f96eca9c41dULL, 0xa3c65b8776f96855ULL, 0x5b34a39f070b5837ULL,
 257     0x4489affce4f31a1eULL, 0x2ffeeb0a48316f40ULL, 0xdc2d9891fe68c022ULL,
 258     0x3659132bb12fea70ULL, 0xaac17d8efa43cab8ULL, 0xc4cb815590989b13ULL,
 259     0x5ee975283d71c93bULL, 0x691548c86c1bd540ULL, 0x7910c41d10a1e6a5ULL,
 260     0x0b5fc64563b3e2a8ULL, 0x047f7684e9fc949dULL, 0xb99181f2d8f685caULL,
 261     0x284600e3f30e38c3ULL
 262   };
 263
 264   uint64_t t[16] = { 0 };
 265   for(size_t i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
 266     for(int b = 0; b < 64; b++)
 267       {
 268         if (JUMP[i] & 1ULL << b)
 269           for(int j = 0; j < 16; j++)
 270             t[j] ^= rs->s[(j + rs->p) & 15];
 271         xorshift1024star (rs);
 272       }
 273   for(int j = 0; j < 16; j++)
 274     rs->s[(j + rs->p) & 15] = t[j];
 275 }
 276
 277
 278 /* Super-simple LCG generator used in getosrandom () if /dev/urandom
 279    doesn't exist.  */
 280
 281 #define M 2147483647 /* 2^31 - 1 (A large prime number) */
 282 #define A 16807      /* Prime root of M, passes statistical tests and produces a full cycle */
 283 #define Q 127773 /* M / A (To avoid overflow on A * seed) */
 284 #define R 2836   /* M % A (To avoid overflow on A * seed) */
 285
 286 __attribute__((unused)) static uint32_t
 287 lcg_parkmiller(uint32_t seed)
 288 {
 289     uint32_t hi = seed / Q;
 290     uint32_t lo = seed % Q;
 291     int32_t test = A * lo - R * hi;
 292     if (test <= 0)
 293         test += M;
 294     return test;
 295 }
 296
 297 #undef M
 298 #undef A
 299 #undef Q
 300 #undef R
 301
 302
 303 /* Get some random bytes from the operating system in order to seed
 304    the PRNG.  */
 305
 306 static int
 307 getosrandom (void *buf, size_t buflen)
 308 {
 309   /* rand_s is available in MinGW-w64 but not plain MinGW.  */
 310 #if defined(__MINGW64_VERSION_MAJOR)
 311   unsigned int* b = buf;
 312   for (size_t i = 0; i < buflen / sizeof (unsigned int); i++)
 313     rand_s (&b[i]);
 314   return buflen;
 315 #else
 316 #ifdef HAVE_GETENTROPY
 317   if (getentropy (buf, buflen) == 0)
 318     return 0;
 319 #endif
 320   int flags = O_RDONLY;
 321 #ifdef O_CLOEXEC
 322   flags |= O_CLOEXEC;
 323 #endif
 324   int fd = open("/dev/urandom", flags);
 325   if (fd != -1)
 326     {
 327       int res = read(fd, buf, buflen);
 328       close (fd);
 329       return res;
 330     }
 331   uint32_t seed = 1234567890;
 332   time_t secs;
 333   long usecs;
 334   if (gf_gettime (&secs, &usecs) == 0)
 335     {
 336       seed ^= secs;
 337       seed ^= usecs;
 338     }
 339 #ifdef HAVE_GETPID
 340   pid_t pid = getpid();
 341   seed ^= pid;
 342 #endif
 343   uint32_t* ub = buf;
 344   for (size_t i = 0; i < buflen / sizeof (uint32_t); i++)
 345     {
 346       ub[i] = seed;
 347       seed = lcg_parkmiller (seed);
 348     }
 349   return buflen;
 350 #endif /* __MINGW64_VERSION_MAJOR  */
 351 }
 352
 353
 354 /* Initialize the random number generator for the current thread,
 355    using the master state and the number of times we must jump.  */
 356
 357 static void
 358 init_rand_state (xorshift1024star_state* rs, const bool locked)
 359 {
 360   if (!locked)
 361     __gthread_mutex_lock (&random_lock);
 362   if (!master_init)
 363     {
 364       getosrandom (master_state, sizeof (master_state));
 365       njumps = 0;
 366       master_init = true;
 367     }
 368   memcpy (&rs->s, master_state, sizeof (master_state));
 369   unsigned n = njumps++;
 370   if (!locked)
 371     __gthread_mutex_unlock (&random_lock);
 372   for (unsigned i = 0; i < n; i++)
 373     jump (rs);
 374   rs->init = true;
 375 }
 376
 377
 378 /*  This function produces a REAL(4) value from the uniform distribution
 379     with range [0,1).  */
 380
 381 void
 382 random_r4 (GFC_REAL_4 *x)
 383 {
 384   xorshift1024star_state* rs = get_rand_state();
 385
 386   if (unlikely (!rs->init))
 387     init_rand_state (rs, false);
 388   uint64_t r = xorshift1024star (rs);
 389   /* Take the higher bits, ensuring that a stream of real(4), real(8),
 390      and real(10) will be identical (except for precision).  */
 391   uint32_t high = (uint32_t) (r >> 32);
 392   rnumber_4 (x, high);
 393 }
 394 iexport(random_r4);
 395
 396 /*  This function produces a REAL(8) value from the uniform distribution
 397     with range [0,1).  */
 398
 399 void
 400 random_r8 (GFC_REAL_8 *x)
 401 {
 402   GFC_UINTEGER_8 r;
 403   xorshift1024star_state* rs = get_rand_state();
 404
 405   if (unlikely (!rs->init))
 406     init_rand_state (rs, false);
 407   r = xorshift1024star (rs);
 408   rnumber_8 (x, r);
 409 }
 410 iexport(random_r8);
 411
 412 #ifdef HAVE_GFC_REAL_10
 413
 414 /*  This function produces a REAL(10) value from the uniform distribution
 415     with range [0,1).  */
 416
 417 void
 418 random_r10 (GFC_REAL_10 *x)
 419 {
 420   GFC_UINTEGER_8 r;
 421   xorshift1024star_state* rs = get_rand_state();
 422
 423   if (unlikely (!rs->init))
 424     init_rand_state (rs, false);
 425   r = xorshift1024star (rs);
 426   rnumber_10 (x, r);
 427 }
 428 iexport(random_r10);
 429
 430 #endif
 431
 432 /*  This function produces a REAL(16) value from the uniform distribution
 433     with range [0,1).  */
 434
 435 #ifdef HAVE_GFC_REAL_16
 436
 437 void
 438 random_r16 (GFC_REAL_16 *x)
 439 {
 440   GFC_UINTEGER_8 r1, r2;
 441   xorshift1024star_state* rs = get_rand_state();
 442
 443   if (unlikely (!rs->init))
 444     init_rand_state (rs, false);
 445   r1 = xorshift1024star (rs);
 446   r2 = xorshift1024star (rs);
 447   rnumber_16 (x, r1, r2);
 448 }
 449 iexport(random_r16);
 450
 451
 452 #endif
 453
 454 /*  This function fills a REAL(4) array with values from the uniform
 455     distribution with range [0,1).  */
 456
 457 void
 458 arandom_r4 (gfc_array_r4 *x)
 459 {
 460   index_type count[GFC_MAX_DIMENSIONS];
 461   index_type extent[GFC_MAX_DIMENSIONS];
 462   index_type stride[GFC_MAX_DIMENSIONS];
 463   index_type stride0;
 464   index_type dim;
 465   GFC_REAL_4 *dest;
 466   xorshift1024star_state* rs = get_rand_state();
 467
 468   dest = x->base_addr;
 469
 470   dim = GFC_DESCRIPTOR_RANK (x);
 471
 472   for (index_type n = 0; n < dim; n++)
 473     {
 474       count[n] = 0;
 475       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 476       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 477       if (extent[n] <= 0)
 478         return;
 479     }
 480
 481   stride0 = stride[0];
 482
 483   if (unlikely (!rs->init))
 484     init_rand_state (rs, false);
 485
 486   while (dest)
 487     {
 488       /* random_r4 (dest);  */
 489       uint64_t r = xorshift1024star (rs);
 490       uint32_t high = (uint32_t) (r >> 32);
 491       rnumber_4 (dest, high);
 492
 493       /* Advance to the next element.  */
 494       dest += stride0;
 495       count[0]++;
 496       /* Advance to the next source element.  */
 497       index_type n = 0;
 498       while (count[n] == extent[n])
 499         {
 500           /* When we get to the end of a dimension, reset it and increment
 501              the next dimension.  */
 502           count[n] = 0;
 503           /* We could precalculate these products, but this is a less
 504              frequently used path so probably not worth it.  */
 505           dest -= stride[n] * extent[n];
 506           n++;
 507           if (n == dim)
 508             {
 509               dest = NULL;
 510               break;
 511             }
 512           else
 513             {
 514               count[n]++;
 515               dest += stride[n];
 516             }
 517         }
 518     }
 519 }
 520
 521 /*  This function fills a REAL(8) array with values from the uniform
 522     distribution with range [0,1).  */
 523
 524 void
 525 arandom_r8 (gfc_array_r8 *x)
 526 {
 527   index_type count[GFC_MAX_DIMENSIONS];
 528   index_type extent[GFC_MAX_DIMENSIONS];
 529   index_type stride[GFC_MAX_DIMENSIONS];
 530   index_type stride0;
 531   index_type dim;
 532   GFC_REAL_8 *dest;
 533   xorshift1024star_state* rs = get_rand_state();
 534
 535   dest = x->base_addr;
 536
 537   dim = GFC_DESCRIPTOR_RANK (x);
 538
 539   for (index_type n = 0; n < dim; n++)
 540     {
 541       count[n] = 0;
 542       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 543       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 544       if (extent[n] <= 0)
 545         return;
 546     }
 547
 548   stride0 = stride[0];
 549
 550   if (unlikely (!rs->init))
 551     init_rand_state (rs, false);
 552
 553   while (dest)
 554     {
 555       /* random_r8 (dest);  */
 556       uint64_t r = xorshift1024star (rs);
 557       rnumber_8 (dest, r);
 558
 559       /* Advance to the next element.  */
 560       dest += stride0;
 561       count[0]++;
 562       /* Advance to the next source element.  */
 563       index_type n = 0;
 564       while (count[n] == extent[n])
 565         {
 566           /* When we get to the end of a dimension, reset it and increment
 567              the next dimension.  */
 568           count[n] = 0;
 569           /* We could precalculate these products, but this is a less
 570              frequently used path so probably not worth it.  */
 571           dest -= stride[n] * extent[n];
 572           n++;
 573           if (n == dim)
 574             {
 575               dest = NULL;
 576               break;
 577             }
 578           else
 579             {
 580               count[n]++;
 581               dest += stride[n];
 582             }
 583         }
 584     }
 585 }
 586
 587 #ifdef HAVE_GFC_REAL_10
 588
 589 /*  This function fills a REAL(10) array with values from the uniform
 590     distribution with range [0,1).  */
 591
 592 void
 593 arandom_r10 (gfc_array_r10 *x)
 594 {
 595   index_type count[GFC_MAX_DIMENSIONS];
 596   index_type extent[GFC_MAX_DIMENSIONS];
 597   index_type stride[GFC_MAX_DIMENSIONS];
 598   index_type stride0;
 599   index_type dim;
 600   GFC_REAL_10 *dest;
 601   xorshift1024star_state* rs = get_rand_state();
 602
 603   dest = x->base_addr;
 604
 605   dim = GFC_DESCRIPTOR_RANK (x);
 606
 607   for (index_type n = 0; n < dim; n++)
 608     {
 609       count[n] = 0;
 610       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 611       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 612       if (extent[n] <= 0)
 613         return;
 614     }
 615
 616   stride0 = stride[0];
 617
 618   if (unlikely (!rs->init))
 619     init_rand_state (rs, false);
 620
 621   while (dest)
 622     {
 623       /* random_r10 (dest);  */
 624       uint64_t r = xorshift1024star (rs);
 625       rnumber_10 (dest, r);
 626
 627       /* Advance to the next element.  */
 628       dest += stride0;
 629       count[0]++;
 630       /* Advance to the next source element.  */
 631       index_type n = 0;
 632       while (count[n] == extent[n])
 633         {
 634           /* When we get to the end of a dimension, reset it and increment
 635              the next dimension.  */
 636           count[n] = 0;
 637           /* We could precalculate these products, but this is a less
 638              frequently used path so probably not worth it.  */
 639           dest -= stride[n] * extent[n];
 640           n++;
 641           if (n == dim)
 642             {
 643               dest = NULL;
 644               break;
 645             }
 646           else
 647             {
 648               count[n]++;
 649               dest += stride[n];
 650             }
 651         }
 652     }
 653 }
 654
 655 #endif
 656
 657 #ifdef HAVE_GFC_REAL_16
 658
 659 /*  This function fills a REAL(16) array with values from the uniform
 660     distribution with range [0,1).  */
 661
 662 void
 663 arandom_r16 (gfc_array_r16 *x)
 664 {
 665   index_type count[GFC_MAX_DIMENSIONS];
 666   index_type extent[GFC_MAX_DIMENSIONS];
 667   index_type stride[GFC_MAX_DIMENSIONS];
 668   index_type stride0;
 669   index_type dim;
 670   GFC_REAL_16 *dest;
 671   xorshift1024star_state* rs = get_rand_state();
 672
 673   dest = x->base_addr;
 674
 675   dim = GFC_DESCRIPTOR_RANK (x);
 676
 677   for (index_type n = 0; n < dim; n++)
 678     {
 679       count[n] = 0;
 680       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 681       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 682       if (extent[n] <= 0)
 683         return;
 684     }
 685
 686   stride0 = stride[0];
 687
 688   if (unlikely (!rs->init))
 689     init_rand_state (rs, false);
 690
 691   while (dest)
 692     {
 693       /* random_r16 (dest);  */
 694       uint64_t r1 = xorshift1024star (rs);
 695       uint64_t r2 = xorshift1024star (rs);
 696       rnumber_16 (dest, r1, r2);
 697
 698       /* Advance to the next element.  */
 699       dest += stride0;
 700       count[0]++;
 701       /* Advance to the next source element.  */
 702       index_type n = 0;
 703       while (count[n] == extent[n])
 704         {
 705           /* When we get to the end of a dimension, reset it and increment
 706              the next dimension.  */
 707           count[n] = 0;
 708           /* We could precalculate these products, but this is a less
 709              frequently used path so probably not worth it.  */
 710           dest -= stride[n] * extent[n];
 711           n++;
 712           if (n == dim)
 713             {
 714               dest = NULL;
 715               break;
 716             }
 717           else
 718             {
 719               count[n]++;
 720               dest += stride[n];
 721             }
 722         }
 723     }
 724 }
 725
 726 #endif
 727
 728
 729 /* Number of elements in master_state array.  */
 730 #define SZU64 (sizeof (master_state) / sizeof (uint64_t))
 731
 732
 733 /* Keys for scrambling the seed in order to avoid poor seeds.  */
 734
 735 static const uint64_t xor_keys[] = {
 736   0xbd0c5b6e50c2df49ULL, 0xd46061cd46e1df38ULL, 0xbb4f4d4ed6103544ULL,
 737   0x114a583d0756ad39ULL, 0x4b5ad8623d0aaab6ULL, 0x3f2ed7afbe0c0f21ULL,
 738   0xdec83fd65f113445ULL, 0x3824f8fbc4f10d24ULL, 0x5d9025af05878911ULL,
 739   0x500bc46b540340e9ULL, 0x8bd53298e0d00530ULL, 0x57886e40a952e06aULL,
 740   0x926e76c88e31cdb6ULL, 0xbd0724dac0a3a5f9ULL, 0xc5c8981b858ab796ULL,
 741   0xbb12ab2694c2b32cULL
 742 };
 743
 744
 745 /* Since a XOR cipher is symmetric, we need only one routine, and we
 746    can use it both for encryption and decryption.  */
 747
 748 static void
 749 scramble_seed (uint64_t *dest, const uint64_t *src)
 750 {
 751   for (size_t i = 0; i < SZU64; i++)
 752     dest[i] = src[i] ^ xor_keys[i];
 753 }
 754
 755
 756 /* random_seed is used to seed the PRNG with either a default
 757    set of seeds or user specified set of seeds.  random_seed
 758    must be called with no argument or exactly one argument.  */
 759
 760 void
 761 random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 762 {
 763   uint64_t seed[SZU64];
 764 #define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_4))
 765
 766   /* Check that we only have one argument present.  */
 767   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 768     runtime_error ("RANDOM_SEED should have at most one argument present.");
 769
 770   if (size != NULL)
 771     *size = SZ + 1;
 772
 773   xorshift1024star_state* rs = get_rand_state();
 774
 775   /* Return the seed to GET data.  */
 776   if (get != NULL)
 777     {
 778       /* If the rank of the array is not 1, abort.  */
 779       if (GFC_DESCRIPTOR_RANK (get) != 1)
 780         runtime_error ("Array rank of GET is not 1.");
 781
 782       /* If the array is too small, abort.  */
 783       if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1)
 784         runtime_error ("Array size of GET is too small.");
 785
 786       if (!rs->init)
 787         init_rand_state (rs, false);
 788
 789       /* Unscramble the seed.  */
 790       scramble_seed (seed, rs->s);
 791
 792       /*  Then copy it back to the user variable.  */
 793       for (size_t i = 0; i < SZ ; i++)
 794         memcpy (&(get->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]),
 795                 (unsigned char*) seed + i * sizeof(GFC_UINTEGER_4),
 796                sizeof(GFC_UINTEGER_4));
 797
 798       /* Finally copy the value of p after the seed.  */
 799       get->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(get, 0)] = rs->p;
 800     }
 801
 802   else
 803     {
 804   __gthread_mutex_lock (&random_lock);
 805
 806   /* From the standard: "If no argument is present, the processor assigns
 807      a processor-dependent value to the seed."  */
 808   if (size == NULL && put == NULL && get == NULL)
 809     {
 810       master_init = false;
 811       init_rand_state (rs, true);
 812     }
 813
 814   if (put != NULL)
 815     {
 816       /* If the rank of the array is not 1, abort.  */
 817       if (GFC_DESCRIPTOR_RANK (put) != 1)
 818         runtime_error ("Array rank of PUT is not 1.");
 819
 820       /* If the array is too small, abort.  */
 821       if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ + 1)
 822         runtime_error ("Array size of PUT is too small.");
 823
 824       /*  We copy the seed given by the user.  */
 825       for (size_t i = 0; i < SZ; i++)
 826         memcpy ((unsigned char*) seed + i * sizeof(GFC_UINTEGER_4),
 827                 &(put->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]),
 828                 sizeof(GFC_UINTEGER_4));
 829
 830       /* We put it after scrambling the bytes, to paper around users who
 831          provide seeds with quality only in the lower or upper part.  */
 832       scramble_seed (master_state, seed);
 833       njumps = 0;
 834       master_init = true;
 835       init_rand_state (rs, true);
 836
 837       rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15;
 838     }
 839
 840   __gthread_mutex_unlock (&random_lock);
 841     }
 842 #undef SZ
 843 }
 844 iexport(random_seed_i4);
 845
 846
 847 void
 848 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
 849 {
 850   uint64_t seed[SZU64];
 851
 852   /* Check that we only have one argument present.  */
 853   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 854     runtime_error ("RANDOM_SEED should have at most one argument present.");
 855
 856 #define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_8))
 857   if (size != NULL)
 858     *size = SZ + 1;
 859
 860   xorshift1024star_state* rs = get_rand_state();
 861
 862   /* Return the seed to GET data.  */
 863   if (get != NULL)
 864     {
 865       /* If the rank of the array is not 1, abort.  */
 866       if (GFC_DESCRIPTOR_RANK (get) != 1)
 867         runtime_error ("Array rank of GET is not 1.");
 868
 869       /* If the array is too small, abort.  */
 870       if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1)
 871         runtime_error ("Array size of GET is too small.");
 872
 873       if (!rs->init)
 874         init_rand_state (rs, false);
 875
 876       /* Unscramble the seed.  */
 877       scramble_seed (seed, rs->s);
 878
 879       /*  This code now should do correct strides.  */
 880       for (size_t i = 0; i < SZ; i++)
 881         memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &seed[i],
 882                 sizeof (GFC_UINTEGER_8));
 883
 884       get->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(get, 0)] = rs->p;
 885     }
 886
 887   else
 888     {
 889   __gthread_mutex_lock (&random_lock);
 890
 891   /* From the standard: "If no argument is present, the processor assigns
 892      a processor-dependent value to the seed."  */
 893   if (size == NULL && put == NULL && get == NULL)
 894     {
 895       master_init = false;
 896       init_rand_state (rs, true);
 897     }
 898
 899   if (put != NULL)
 900     {
 901       /* If the rank of the array is not 1, abort.  */
 902       if (GFC_DESCRIPTOR_RANK (put) != 1)
 903         runtime_error ("Array rank of PUT is not 1.");
 904
 905       /* If the array is too small, abort.  */
 906       if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ + 1)
 907         runtime_error ("Array size of PUT is too small.");
 908
 909       /*  This code now should do correct strides.  */
 910       for (size_t i = 0; i < SZ; i++)
 911         memcpy (&seed[i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]),
 912                 sizeof (GFC_UINTEGER_8));
 913
 914       scramble_seed (master_state, seed);
 915       njumps = 0;
 916       master_init = true;
 917       init_rand_state (rs, true);
 918       rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15;
 919      }
 920
 921
 922   __gthread_mutex_unlock (&random_lock);
 923     }
 924 }
 925 iexport(random_seed_i8);
 926
 927
 928 #if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS
 929 static void __attribute__((constructor))
 930 constructor_random (void)
 931 {
 932 #ifndef __GTHREAD_MUTEX_INIT
 933   __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
 934 #endif
 935   if (__gthread_active_p ())
 936     __gthread_key_create (&rand_state_key, &free);
 937 }
 938 #endif
 939
 940 #ifdef __GTHREADS
 941 static void __attribute__((destructor))
 942 destructor_random (void)
 943 {
 944   if (__gthread_active_p ())
 945     __gthread_key_delete (rand_state_key);
 946 }
 947 #endif