libgfortran/intrinsics/random.c

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