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 #ifdef HAVE_GETENTROPY
 314   if (getentropy (buf, buflen) == 0)
 315     return 0;
 316 #endif
 317   int flags = O_RDONLY;
 318 #ifdef O_CLOEXEC
 319   flags |= O_CLOEXEC;
 320 #endif
 321   int fd = open("/dev/urandom", flags);
 322   if (fd != -1)
 323     {
 324       int res = read(fd, buf, buflen);
 325       close (fd);
 326       return res;
 327     }
 328   uint32_t seed = 1234567890;
 329   time_t secs;
 330   long usecs;
 331   if (gf_gettime (&secs, &usecs) == 0)
 332     {
 333       seed ^= secs;
 334       seed ^= usecs;
 335     }
 336 #ifdef HAVE_GETPID
 337   pid_t pid = getpid();
 338   seed ^= pid;
 339 #endif
 340   uint32_t* ub = buf;
 341   for (size_t i = 0; i < buflen / sizeof (uint32_t); i++)
 342     {
 343       ub[i] = seed;
 344       seed = lcg_parkmiller (seed);
 345     }
 346   return buflen;
 347 #endif /* __MINGW64_VERSION_MAJOR  */
 348 }
 349
 350
 351 /* Initialize the random number generator for the current thread,
 352    using the master state and the number of times we must jump.  */
 353
 354 static void
 355 init_rand_state (xorshift1024star_state* rs, const bool locked)
 356 {
 357   if (!locked)
 358     __gthread_mutex_lock (&random_lock);
 359   if (!master_init)
 360     {
 361       getosrandom (master_state, sizeof (master_state));
 362       njumps = 0;
 363       master_init = true;
 364     }
 365   memcpy (&rs->s, master_state, sizeof (master_state));
 366   unsigned n = njumps++;
 367   if (!locked)
 368     __gthread_mutex_unlock (&random_lock);
 369   for (unsigned i = 0; i < n; i++)
 370     jump (rs);
 371   rs->init = true;
 372 }
 373
 374
 375 /*  This function produces a REAL(4) value from the uniform distribution
 376     with range [0,1).  */
 377
 378 void
 379 random_r4 (GFC_REAL_4 *x)
 380 {
 381   xorshift1024star_state* rs = get_rand_state();
 382
 383   if (unlikely (!rs->init))
 384     init_rand_state (rs, false);
 385   uint64_t r = xorshift1024star (rs);
 386   /* Take the higher bits, ensuring that a stream of real(4), real(8),
 387      and real(10) will be identical (except for precision).  */
 388   uint32_t high = (uint32_t) (r >> 32);
 389   rnumber_4 (x, high);
 390 }
 391 iexport(random_r4);
 392
 393 /*  This function produces a REAL(8) value from the uniform distribution
 394     with range [0,1).  */
 395
 396 void
 397 random_r8 (GFC_REAL_8 *x)
 398 {
 399   GFC_UINTEGER_8 r;
 400   xorshift1024star_state* rs = get_rand_state();
 401
 402   if (unlikely (!rs->init))
 403     init_rand_state (rs, false);
 404   r = xorshift1024star (rs);
 405   rnumber_8 (x, r);
 406 }
 407 iexport(random_r8);
 408
 409 #ifdef HAVE_GFC_REAL_10
 410
 411 /*  This function produces a REAL(10) value from the uniform distribution
 412     with range [0,1).  */
 413
 414 void
 415 random_r10 (GFC_REAL_10 *x)
 416 {
 417   GFC_UINTEGER_8 r;
 418   xorshift1024star_state* rs = get_rand_state();
 419
 420   if (unlikely (!rs->init))
 421     init_rand_state (rs, false);
 422   r = xorshift1024star (rs);
 423   rnumber_10 (x, r);
 424 }
 425 iexport(random_r10);
 426
 427 #endif
 428
 429 /*  This function produces a REAL(16) value from the uniform distribution
 430     with range [0,1).  */
 431
 432 #ifdef HAVE_GFC_REAL_16
 433
 434 void
 435 random_r16 (GFC_REAL_16 *x)
 436 {
 437   GFC_UINTEGER_8 r1, r2;
 438   xorshift1024star_state* rs = get_rand_state();
 439
 440   if (unlikely (!rs->init))
 441     init_rand_state (rs, false);
 442   r1 = xorshift1024star (rs);
 443   r2 = xorshift1024star (rs);
 444   rnumber_16 (x, r1, r2);
 445 }
 446 iexport(random_r16);
 447
 448
 449 #endif
 450
 451 /*  This function fills a REAL(4) array with values from the uniform
 452     distribution with range [0,1).  */
 453
 454 void
 455 arandom_r4 (gfc_array_r4 *x)
 456 {
 457   index_type count[GFC_MAX_DIMENSIONS];
 458   index_type extent[GFC_MAX_DIMENSIONS];
 459   index_type stride[GFC_MAX_DIMENSIONS];
 460   index_type stride0;
 461   index_type dim;
 462   GFC_REAL_4 *dest;
 463   xorshift1024star_state* rs = get_rand_state();
 464
 465   dest = x->base_addr;
 466
 467   dim = GFC_DESCRIPTOR_RANK (x);
 468
 469   for (index_type n = 0; n < dim; n++)
 470     {
 471       count[n] = 0;
 472       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 473       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 474       if (extent[n] <= 0)
 475         return;
 476     }
 477
 478   stride0 = stride[0];
 479
 480   if (unlikely (!rs->init))
 481     init_rand_state (rs, false);
 482
 483   while (dest)
 484     {
 485       /* random_r4 (dest);  */
 486       uint64_t r = xorshift1024star (rs);
 487       uint32_t high = (uint32_t) (r >> 32);
 488       rnumber_4 (dest, high);
 489
 490       /* Advance to the next element.  */
 491       dest += stride0;
 492       count[0]++;
 493       /* Advance to the next source element.  */
 494       index_type n = 0;
 495       while (count[n] == extent[n])
 496         {
 497           /* When we get to the end of a dimension, reset it and increment
 498              the next dimension.  */
 499           count[n] = 0;
 500           /* We could precalculate these products, but this is a less
 501              frequently used path so probably not worth it.  */
 502           dest -= stride[n] * extent[n];
 503           n++;
 504           if (n == dim)
 505             {
 506               dest = NULL;
 507               break;
 508             }
 509           else
 510             {
 511               count[n]++;
 512               dest += stride[n];
 513             }
 514         }
 515     }
 516 }
 517
 518 /*  This function fills a REAL(8) array with values from the uniform
 519     distribution with range [0,1).  */
 520
 521 void
 522 arandom_r8 (gfc_array_r8 *x)
 523 {
 524   index_type count[GFC_MAX_DIMENSIONS];
 525   index_type extent[GFC_MAX_DIMENSIONS];
 526   index_type stride[GFC_MAX_DIMENSIONS];
 527   index_type stride0;
 528   index_type dim;
 529   GFC_REAL_8 *dest;
 530   xorshift1024star_state* rs = get_rand_state();
 531
 532   dest = x->base_addr;
 533
 534   dim = GFC_DESCRIPTOR_RANK (x);
 535
 536   for (index_type n = 0; n < dim; n++)
 537     {
 538       count[n] = 0;
 539       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 540       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 541       if (extent[n] <= 0)
 542         return;
 543     }
 544
 545   stride0 = stride[0];
 546
 547   if (unlikely (!rs->init))
 548     init_rand_state (rs, false);
 549
 550   while (dest)
 551     {
 552       /* random_r8 (dest);  */
 553       uint64_t r = xorshift1024star (rs);
 554       rnumber_8 (dest, r);
 555
 556       /* Advance to the next element.  */
 557       dest += stride0;
 558       count[0]++;
 559       /* Advance to the next source element.  */
 560       index_type n = 0;
 561       while (count[n] == extent[n])
 562         {
 563           /* When we get to the end of a dimension, reset it and increment
 564              the next dimension.  */
 565           count[n] = 0;
 566           /* We could precalculate these products, but this is a less
 567              frequently used path so probably not worth it.  */
 568           dest -= stride[n] * extent[n];
 569           n++;
 570           if (n == dim)
 571             {
 572               dest = NULL;
 573               break;
 574             }
 575           else
 576             {
 577               count[n]++;
 578               dest += stride[n];
 579             }
 580         }
 581     }
 582 }
 583
 584 #ifdef HAVE_GFC_REAL_10
 585
 586 /*  This function fills a REAL(10) array with values from the uniform
 587     distribution with range [0,1).  */
 588
 589 void
 590 arandom_r10 (gfc_array_r10 *x)
 591 {
 592   index_type count[GFC_MAX_DIMENSIONS];
 593   index_type extent[GFC_MAX_DIMENSIONS];
 594   index_type stride[GFC_MAX_DIMENSIONS];
 595   index_type stride0;
 596   index_type dim;
 597   GFC_REAL_10 *dest;
 598   xorshift1024star_state* rs = get_rand_state();
 599
 600   dest = x->base_addr;
 601
 602   dim = GFC_DESCRIPTOR_RANK (x);
 603
 604   for (index_type n = 0; n < dim; n++)
 605     {
 606       count[n] = 0;
 607       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 608       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 609       if (extent[n] <= 0)
 610         return;
 611     }
 612
 613   stride0 = stride[0];
 614
 615   if (unlikely (!rs->init))
 616     init_rand_state (rs, false);
 617
 618   while (dest)
 619     {
 620       /* random_r10 (dest);  */
 621       uint64_t r = xorshift1024star (rs);
 622       rnumber_10 (dest, r);
 623
 624       /* Advance to the next element.  */
 625       dest += stride0;
 626       count[0]++;
 627       /* Advance to the next source element.  */
 628       index_type n = 0;
 629       while (count[n] == extent[n])
 630         {
 631           /* When we get to the end of a dimension, reset it and increment
 632              the next dimension.  */
 633           count[n] = 0;
 634           /* We could precalculate these products, but this is a less
 635              frequently used path so probably not worth it.  */
 636           dest -= stride[n] * extent[n];
 637           n++;
 638           if (n == dim)
 639             {
 640               dest = NULL;
 641               break;
 642             }
 643           else
 644             {
 645               count[n]++;
 646               dest += stride[n];
 647             }
 648         }
 649     }
 650 }
 651
 652 #endif
 653
 654 #ifdef HAVE_GFC_REAL_16
 655
 656 /*  This function fills a REAL(16) array with values from the uniform
 657     distribution with range [0,1).  */
 658
 659 void
 660 arandom_r16 (gfc_array_r16 *x)
 661 {
 662   index_type count[GFC_MAX_DIMENSIONS];
 663   index_type extent[GFC_MAX_DIMENSIONS];
 664   index_type stride[GFC_MAX_DIMENSIONS];
 665   index_type stride0;
 666   index_type dim;
 667   GFC_REAL_16 *dest;
 668   xorshift1024star_state* rs = get_rand_state();
 669
 670   dest = x->base_addr;
 671
 672   dim = GFC_DESCRIPTOR_RANK (x);
 673
 674   for (index_type n = 0; n < dim; n++)
 675     {
 676       count[n] = 0;
 677       stride[n] = GFC_DESCRIPTOR_STRIDE(x,n);
 678       extent[n] = GFC_DESCRIPTOR_EXTENT(x,n);
 679       if (extent[n] <= 0)
 680         return;
 681     }
 682
 683   stride0 = stride[0];
 684
 685   if (unlikely (!rs->init))
 686     init_rand_state (rs, false);
 687
 688   while (dest)
 689     {
 690       /* random_r16 (dest);  */
 691       uint64_t r1 = xorshift1024star (rs);
 692       uint64_t r2 = xorshift1024star (rs);
 693       rnumber_16 (dest, r1, r2);
 694
 695       /* Advance to the next element.  */
 696       dest += stride0;
 697       count[0]++;
 698       /* Advance to the next source element.  */
 699       index_type n = 0;
 700       while (count[n] == extent[n])
 701         {
 702           /* When we get to the end of a dimension, reset it and increment
 703              the next dimension.  */
 704           count[n] = 0;
 705           /* We could precalculate these products, but this is a less
 706              frequently used path so probably not worth it.  */
 707           dest -= stride[n] * extent[n];
 708           n++;
 709           if (n == dim)
 710             {
 711               dest = NULL;
 712               break;
 713             }
 714           else
 715             {
 716               count[n]++;
 717               dest += stride[n];
 718             }
 719         }
 720     }
 721 }
 722
 723 #endif
 724
 725
 726 /* Number of elements in master_state array.  */
 727 #define SZU64 (sizeof (master_state) / sizeof (uint64_t))
 728
 729
 730 /* Keys for scrambling the seed in order to avoid poor seeds.  */
 731
 732 static const uint64_t xor_keys[] = {
 733   0xbd0c5b6e50c2df49ULL, 0xd46061cd46e1df38ULL, 0xbb4f4d4ed6103544ULL,
 734   0x114a583d0756ad39ULL, 0x4b5ad8623d0aaab6ULL, 0x3f2ed7afbe0c0f21ULL,
 735   0xdec83fd65f113445ULL, 0x3824f8fbc4f10d24ULL, 0x5d9025af05878911ULL,
 736   0x500bc46b540340e9ULL, 0x8bd53298e0d00530ULL, 0x57886e40a952e06aULL,
 737   0x926e76c88e31cdb6ULL, 0xbd0724dac0a3a5f9ULL, 0xc5c8981b858ab796ULL,
 738   0xbb12ab2694c2b32cULL
 739 };
 740
 741
 742 /* Since a XOR cipher is symmetric, we need only one routine, and we
 743    can use it both for encryption and decryption.  */
 744
 745 static void
 746 scramble_seed (uint64_t *dest, const uint64_t *src)
 747 {
 748   for (size_t i = 0; i < SZU64; i++)
 749     dest[i] = src[i] ^ xor_keys[i];
 750 }
 751
 752
 753 /* random_seed is used to seed the PRNG with either a default
 754    set of seeds or user specified set of seeds.  random_seed
 755    must be called with no argument or exactly one argument.  */
 756
 757 void
 758 random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
 759 {
 760   uint64_t seed[SZU64];
 761 #define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_4))
 762
 763   /* Check that we only have one argument present.  */
 764   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 765     runtime_error ("RANDOM_SEED should have at most one argument present.");
 766
 767   if (size != NULL)
 768     *size = SZ + 1;
 769
 770   xorshift1024star_state* rs = get_rand_state();
 771
 772   /* Return the seed to GET data.  */
 773   if (get != NULL)
 774     {
 775       /* If the rank of the array is not 1, abort.  */
 776       if (GFC_DESCRIPTOR_RANK (get) != 1)
 777         runtime_error ("Array rank of GET is not 1.");
 778
 779       /* If the array is too small, abort.  */
 780       if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1)
 781         runtime_error ("Array size of GET is too small.");
 782
 783       if (!rs->init)
 784         init_rand_state (rs, false);
 785
 786       /* Unscramble the seed.  */
 787       scramble_seed (seed, rs->s);
 788
 789       /*  Then copy it back to the user variable.  */
 790       for (size_t i = 0; i < SZ ; i++)
 791         memcpy (&(get->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]),
 792                 (unsigned char*) seed + i * sizeof(GFC_UINTEGER_4),
 793                sizeof(GFC_UINTEGER_4));
 794
 795       /* Finally copy the value of p after the seed.  */
 796       get->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(get, 0)] = rs->p;
 797     }
 798
 799   else
 800     {
 801   __gthread_mutex_lock (&random_lock);
 802
 803   /* From the standard: "If no argument is present, the processor assigns
 804      a processor-dependent value to the seed."  */
 805   if (size == NULL && put == NULL && get == NULL)
 806     {
 807       master_init = false;
 808       init_rand_state (rs, true);
 809     }
 810
 811   if (put != NULL)
 812     {
 813       /* If the rank of the array is not 1, abort.  */
 814       if (GFC_DESCRIPTOR_RANK (put) != 1)
 815         runtime_error ("Array rank of PUT is not 1.");
 816
 817       /* If the array is too small, abort.  */
 818       if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ + 1)
 819         runtime_error ("Array size of PUT is too small.");
 820
 821       /*  We copy the seed given by the user.  */
 822       for (size_t i = 0; i < SZ; i++)
 823         memcpy ((unsigned char*) seed + i * sizeof(GFC_UINTEGER_4),
 824                 &(put->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]),
 825                 sizeof(GFC_UINTEGER_4));
 826
 827       /* We put it after scrambling the bytes, to paper around users who
 828          provide seeds with quality only in the lower or upper part.  */
 829       scramble_seed (master_state, seed);
 830       njumps = 0;
 831       master_init = true;
 832       init_rand_state (rs, true);
 833
 834       rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15;
 835     }
 836
 837   __gthread_mutex_unlock (&random_lock);
 838     }
 839 #undef SZ
 840 }
 841 iexport(random_seed_i4);
 842
 843
 844 void
 845 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get)
 846 {
 847   uint64_t seed[SZU64];
 848
 849   /* Check that we only have one argument present.  */
 850   if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1)
 851     runtime_error ("RANDOM_SEED should have at most one argument present.");
 852
 853 #define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_8))
 854   if (size != NULL)
 855     *size = SZ + 1;
 856
 857   xorshift1024star_state* rs = get_rand_state();
 858
 859   /* Return the seed to GET data.  */
 860   if (get != NULL)
 861     {
 862       /* If the rank of the array is not 1, abort.  */
 863       if (GFC_DESCRIPTOR_RANK (get) != 1)
 864         runtime_error ("Array rank of GET is not 1.");
 865
 866       /* If the array is too small, abort.  */
 867       if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ + 1)
 868         runtime_error ("Array size of GET is too small.");
 869
 870       if (!rs->init)
 871         init_rand_state (rs, false);
 872
 873       /* Unscramble the seed.  */
 874       scramble_seed (seed, rs->s);
 875
 876       /*  This code now should do correct strides.  */
 877       for (size_t i = 0; i < SZ; i++)
 878         memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &seed[i],
 879                 sizeof (GFC_UINTEGER_8));
 880
 881       get->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(get, 0)] = rs->p;
 882     }
 883
 884   else
 885     {
 886   __gthread_mutex_lock (&random_lock);
 887
 888   /* From the standard: "If no argument is present, the processor assigns
 889      a processor-dependent value to the seed."  */
 890   if (size == NULL && put == NULL && get == NULL)
 891     {
 892       master_init = false;
 893       init_rand_state (rs, true);
 894     }
 895
 896   if (put != NULL)
 897     {
 898       /* If the rank of the array is not 1, abort.  */
 899       if (GFC_DESCRIPTOR_RANK (put) != 1)
 900         runtime_error ("Array rank of PUT is not 1.");
 901
 902       /* If the array is too small, abort.  */
 903       if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ + 1)
 904         runtime_error ("Array size of PUT is too small.");
 905
 906       /*  This code now should do correct strides.  */
 907       for (size_t i = 0; i < SZ; i++)
 908         memcpy (&seed[i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]),
 909                 sizeof (GFC_UINTEGER_8));
 910
 911       scramble_seed (master_state, seed);
 912       njumps = 0;
 913       master_init = true;
 914       init_rand_state (rs, true);
 915       rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15;
 916      }
 917
 918
 919   __gthread_mutex_unlock (&random_lock);
 920     }
 921 }
 922 iexport(random_seed_i8);
 923
 924
 925 #if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS
 926 static void __attribute__((constructor))
 927 constructor_random (void)
 928 {
 929 #ifndef __GTHREAD_MUTEX_INIT
 930   __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock);
 931 #endif
 932   if (__gthread_active_p ())
 933     __gthread_key_create (&rand_state_key, &free);
 934 }
 935 #endif
 936
 937 #ifdef __GTHREADS
 938 static void __attribute__((destructor))
 939 destructor_random (void)
 940 {
 941   if (__gthread_active_p ())
 942     __gthread_key_delete (rand_state_key);
 943 }
 944 #endif