libgfortran/intrinsics/random.c

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