libgfortran/intrinsics/random.c

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