libgfortran/intrinsics/random.c

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