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