1 /* Implementation of the RANDOM intrinsics 2 Copyright (C) 2002-2019 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 xorshift1024* 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**1024 - 1 176 177 A description can be found at 178 179 http://vigna.di.unimi.it/ftp/papers/xorshift.pdf 180 181 or 182 183 http://arxiv.org/abs/1402.6246 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 int p; 193 uint64_t s[16]; 194 } 195 xorshift1024star_state; 196 197 198 /* master_init, njumps, and master_state are the only variables 199 protected by random_lock. */ 200 static bool master_init; 201 static unsigned njumps; /* How many times we have jumped. */ 202 static uint64_t master_state[] = { 203 0xad63fa1ed3b55f36ULL, 0xd94473e78978b497ULL, 0xbc60592a98172477ULL, 204 0xa3de7c6e81265301ULL, 0x586640c5e785af27ULL, 0x7a2a3f63b67ce5eaULL, 205 0x9fde969f922d9b82ULL, 0xe6fe34379b3f3822ULL, 0x6c277eac3e99b6c2ULL, 206 0x9197290ab0d3f069ULL, 0xdb227302f6c25576ULL, 0xee0209aee527fae9ULL, 207 0x675666a793cd05b9ULL, 0xd048c99fbc70c20fULL, 0x775f8c3dba385ef5ULL, 208 0x625288bc262faf33ULL 209 }; 210 211 212 static __gthread_key_t rand_state_key; 213 214 static xorshift1024star_state* 215 get_rand_state (void) 216 { 217 /* For single threaded apps. */ 218 static xorshift1024star_state rand_state; 219 220 if (__gthread_active_p ()) 221 { 222 void* p = __gthread_getspecific (rand_state_key); 223 if (!p) 224 { 225 p = xcalloc (1, sizeof (xorshift1024star_state)); 226 __gthread_setspecific (rand_state_key, p); 227 } 228 return p; 229 } 230 else 231 return &rand_state; 232 } 233 234 235 static uint64_t 236 xorshift1024star (xorshift1024star_state* rs) 237 { 238 int p = rs->p; 239 const uint64_t s0 = rs->s[p]; 240 uint64_t s1 = rs->s[p = (p + 1) & 15]; 241 s1 ^= s1 << 31; 242 rs->s[p] = s1 ^ s0 ^ (s1 >> 11) ^ (s0 >> 30); 243 rs->p = p; 244 return rs->s[p] * UINT64_C(1181783497276652981); 245 } 246 247 248 /* This is the jump function for the generator. It is equivalent to 249 2^512 calls to xorshift1024star(); it can be used to generate 2^512 250 non-overlapping subsequences for parallel computations. */ 251 252 static void 253 jump (xorshift1024star_state* rs) 254 { 255 static const uint64_t JUMP[] = { 256 0x84242f96eca9c41dULL, 0xa3c65b8776f96855ULL, 0x5b34a39f070b5837ULL, 257 0x4489affce4f31a1eULL, 0x2ffeeb0a48316f40ULL, 0xdc2d9891fe68c022ULL, 258 0x3659132bb12fea70ULL, 0xaac17d8efa43cab8ULL, 0xc4cb815590989b13ULL, 259 0x5ee975283d71c93bULL, 0x691548c86c1bd540ULL, 0x7910c41d10a1e6a5ULL, 260 0x0b5fc64563b3e2a8ULL, 0x047f7684e9fc949dULL, 0xb99181f2d8f685caULL, 261 0x284600e3f30e38c3ULL 262 }; 263 264 uint64_t t[16] = { 0 }; 265 for(size_t i = 0; i < sizeof JUMP / sizeof *JUMP; i++) 266 for(int b = 0; b < 64; b++) 267 { 268 if (JUMP[i] & 1ULL << b) 269 for(int j = 0; j < 16; j++) 270 t[j] ^= rs->s[(j + rs->p) & 15]; 271 xorshift1024star (rs); 272 } 273 for(int j = 0; j < 16; j++) 274 rs->s[(j + rs->p) & 15] = t[j]; 275 } 276 277 278 /* Splitmix64 recommended by xorshift author for initializing. After 279 getting one uint64_t value from the OS, this is used to fill in the 280 rest of the state. */ 281 282 static uint64_t 283 splitmix64 (uint64_t x) 284 { 285 uint64_t z = (x += 0x9e3779b97f4a7c15); 286 z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9; 287 z = (z ^ (z >> 27)) * 0x94d049bb133111eb; 288 return z ^ (z >> 31); 289 } 290 291 292 /* Get some random bytes from the operating system in order to seed 293 the PRNG. */ 294 295 static int 296 getosrandom (void *buf, size_t buflen) 297 { 298 /* rand_s is available in MinGW-w64 but not plain MinGW. */ 299 #if defined(__MINGW64_VERSION_MAJOR) 300 unsigned int* b = buf; 301 for (size_t i = 0; i < buflen / sizeof (unsigned int); i++) 302 rand_s (&b[i]); 303 return buflen; 304 #else 305 #ifdef HAVE_GETENTROPY 306 if (getentropy (buf, buflen) == 0) 307 return buflen; 308 #endif 309 int flags = O_RDONLY; 310 #ifdef O_CLOEXEC 311 flags |= O_CLOEXEC; 312 #endif 313 int fd = open("/dev/urandom", flags); 314 if (fd != -1) 315 { 316 int res = read(fd, buf, buflen); 317 close (fd); 318 return res; 319 } 320 uint64_t seed = 0x047f7684e9fc949dULL; 321 time_t secs; 322 long usecs; 323 if (gf_gettime (&secs, &usecs) == 0) 324 { 325 seed ^= secs; 326 seed ^= usecs; 327 } 328 #ifdef HAVE_GETPID 329 pid_t pid = getpid(); 330 seed ^= pid; 331 #endif 332 size_t size = buflen < sizeof (uint64_t) ? buflen : sizeof (uint64_t); 333 memcpy (buf, &seed, size); 334 return size; 335 #endif /* __MINGW64_VERSION_MAJOR */ 336 } 337 338 339 /* Initialize the random number generator for the current thread, 340 using the master state and the number of times we must jump. */ 341 342 static void 343 init_rand_state (xorshift1024star_state* rs, const bool locked) 344 { 345 if (!locked) 346 __gthread_mutex_lock (&random_lock); 347 if (!master_init) 348 { 349 uint64_t os_seed; 350 getosrandom (&os_seed, sizeof (os_seed)); 351 for (uint64_t i = 0; i < sizeof (master_state) / sizeof (uint64_t); i++) 352 { 353 os_seed = splitmix64 (os_seed); 354 master_state[i] = os_seed; 355 } 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 459 dest = x->base_addr; 460 461 dim = GFC_DESCRIPTOR_RANK (x); 462 463 for (index_type n = 0; n < dim; n++) 464 { 465 count[n] = 0; 466 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n); 467 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n); 468 if (extent[n] <= 0) 469 return; 470 } 471 472 stride0 = stride[0]; 473 474 if (unlikely (!rs->init)) 475 init_rand_state (rs, false); 476 477 while (dest) 478 { 479 /* random_r4 (dest); */ 480 uint64_t r = xorshift1024star (rs); 481 uint32_t high = (uint32_t) (r >> 32); 482 rnumber_4 (dest, high); 483 484 /* Advance to the next element. */ 485 dest += stride0; 486 count[0]++; 487 /* Advance to the next source element. */ 488 index_type n = 0; 489 while (count[n] == extent[n]) 490 { 491 /* When we get to the end of a dimension, reset it and increment 492 the next dimension. */ 493 count[n] = 0; 494 /* We could precalculate these products, but this is a less 495 frequently used path so probably not worth it. */ 496 dest -= stride[n] * extent[n]; 497 n++; 498 if (n == dim) 499 { 500 dest = NULL; 501 break; 502 } 503 else 504 { 505 count[n]++; 506 dest += stride[n]; 507 } 508 } 509 } 510 } 511 512 /* This function fills a REAL(8) array with values from the uniform 513 distribution with range [0,1). */ 514 515 void 516 arandom_r8 (gfc_array_r8 *x) 517 { 518 index_type count[GFC_MAX_DIMENSIONS]; 519 index_type extent[GFC_MAX_DIMENSIONS]; 520 index_type stride[GFC_MAX_DIMENSIONS]; 521 index_type stride0; 522 index_type dim; 523 GFC_REAL_8 *dest; 524 xorshift1024star_state* rs = get_rand_state(); 525 526 dest = x->base_addr; 527 528 dim = GFC_DESCRIPTOR_RANK (x); 529 530 for (index_type n = 0; n < dim; n++) 531 { 532 count[n] = 0; 533 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n); 534 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n); 535 if (extent[n] <= 0) 536 return; 537 } 538 539 stride0 = stride[0]; 540 541 if (unlikely (!rs->init)) 542 init_rand_state (rs, false); 543 544 while (dest) 545 { 546 /* random_r8 (dest); */ 547 uint64_t r = xorshift1024star (rs); 548 rnumber_8 (dest, r); 549 550 /* Advance to the next element. */ 551 dest += stride0; 552 count[0]++; 553 /* Advance to the next source element. */ 554 index_type n = 0; 555 while (count[n] == extent[n]) 556 { 557 /* When we get to the end of a dimension, reset it and increment 558 the next dimension. */ 559 count[n] = 0; 560 /* We could precalculate these products, but this is a less 561 frequently used path so probably not worth it. */ 562 dest -= stride[n] * extent[n]; 563 n++; 564 if (n == dim) 565 { 566 dest = NULL; 567 break; 568 } 569 else 570 { 571 count[n]++; 572 dest += stride[n]; 573 } 574 } 575 } 576 } 577 578 #ifdef HAVE_GFC_REAL_10 579 580 /* This function fills a REAL(10) array with values from the uniform 581 distribution with range [0,1). */ 582 583 void 584 arandom_r10 (gfc_array_r10 *x) 585 { 586 index_type count[GFC_MAX_DIMENSIONS]; 587 index_type extent[GFC_MAX_DIMENSIONS]; 588 index_type stride[GFC_MAX_DIMENSIONS]; 589 index_type stride0; 590 index_type dim; 591 GFC_REAL_10 *dest; 592 xorshift1024star_state* rs = get_rand_state(); 593 594 dest = x->base_addr; 595 596 dim = GFC_DESCRIPTOR_RANK (x); 597 598 for (index_type n = 0; n < dim; n++) 599 { 600 count[n] = 0; 601 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n); 602 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n); 603 if (extent[n] <= 0) 604 return; 605 } 606 607 stride0 = stride[0]; 608 609 if (unlikely (!rs->init)) 610 init_rand_state (rs, false); 611 612 while (dest) 613 { 614 /* random_r10 (dest); */ 615 uint64_t r = xorshift1024star (rs); 616 rnumber_10 (dest, r); 617 618 /* Advance to the next element. */ 619 dest += stride0; 620 count[0]++; 621 /* Advance to the next source element. */ 622 index_type n = 0; 623 while (count[n] == extent[n]) 624 { 625 /* When we get to the end of a dimension, reset it and increment 626 the next dimension. */ 627 count[n] = 0; 628 /* We could precalculate these products, but this is a less 629 frequently used path so probably not worth it. */ 630 dest -= stride[n] * extent[n]; 631 n++; 632 if (n == dim) 633 { 634 dest = NULL; 635 break; 636 } 637 else 638 { 639 count[n]++; 640 dest += stride[n]; 641 } 642 } 643 } 644 } 645 646 #endif 647 648 #ifdef HAVE_GFC_REAL_16 649 650 /* This function fills a REAL(16) array with values from the uniform 651 distribution with range [0,1). */ 652 653 void 654 arandom_r16 (gfc_array_r16 *x) 655 { 656 index_type count[GFC_MAX_DIMENSIONS]; 657 index_type extent[GFC_MAX_DIMENSIONS]; 658 index_type stride[GFC_MAX_DIMENSIONS]; 659 index_type stride0; 660 index_type dim; 661 GFC_REAL_16 *dest; 662 xorshift1024star_state* rs = get_rand_state(); 663 664 dest = x->base_addr; 665 666 dim = GFC_DESCRIPTOR_RANK (x); 667 668 for (index_type n = 0; n < dim; n++) 669 { 670 count[n] = 0; 671 stride[n] = GFC_DESCRIPTOR_STRIDE(x,n); 672 extent[n] = GFC_DESCRIPTOR_EXTENT(x,n); 673 if (extent[n] <= 0) 674 return; 675 } 676 677 stride0 = stride[0]; 678 679 if (unlikely (!rs->init)) 680 init_rand_state (rs, false); 681 682 while (dest) 683 { 684 /* random_r16 (dest); */ 685 uint64_t r1 = xorshift1024star (rs); 686 uint64_t r2 = xorshift1024star (rs); 687 rnumber_16 (dest, r1, r2); 688 689 /* Advance to the next element. */ 690 dest += stride0; 691 count[0]++; 692 /* Advance to the next source element. */ 693 index_type n = 0; 694 while (count[n] == extent[n]) 695 { 696 /* When we get to the end of a dimension, reset it and increment 697 the next dimension. */ 698 count[n] = 0; 699 /* We could precalculate these products, but this is a less 700 frequently used path so probably not worth it. */ 701 dest -= stride[n] * extent[n]; 702 n++; 703 if (n == dim) 704 { 705 dest = NULL; 706 break; 707 } 708 else 709 { 710 count[n]++; 711 dest += stride[n]; 712 } 713 } 714 } 715 } 716 717 #endif 718 719 720 /* Number of elements in master_state array. */ 721 #define SZU64 (sizeof (master_state) / sizeof (uint64_t)) 722 723 724 /* Keys for scrambling the seed in order to avoid poor seeds. */ 725 726 static const uint64_t xor_keys[] = { 727 0xbd0c5b6e50c2df49ULL, 0xd46061cd46e1df38ULL, 0xbb4f4d4ed6103544ULL, 728 0x114a583d0756ad39ULL, 0x4b5ad8623d0aaab6ULL, 0x3f2ed7afbe0c0f21ULL, 729 0xdec83fd65f113445ULL, 0x3824f8fbc4f10d24ULL, 0x5d9025af05878911ULL, 730 0x500bc46b540340e9ULL, 0x8bd53298e0d00530ULL, 0x57886e40a952e06aULL, 731 0x926e76c88e31cdb6ULL, 0xbd0724dac0a3a5f9ULL, 0xc5c8981b858ab796ULL, 732 0xbb12ab2694c2b32cULL 733 }; 734 735 736 /* Since a XOR cipher is symmetric, we need only one routine, and we 737 can use it both for encryption and decryption. */ 738 739 static void 740 scramble_seed (uint64_t *dest, const uint64_t *src) 741 { 742 for (size_t i = 0; i < SZU64; i++) 743 dest[i] = src[i] ^ xor_keys[i]; 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 uint64_t seed[SZU64]; 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 scramble_seed (seed, rs->s); 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 (unsigned char*) 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 ((unsigned char*) 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 (master_state, 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 __gthread_mutex_unlock (&random_lock); 832 } 833 #undef SZ 834 } 835 iexport(random_seed_i4); 836 837 838 void 839 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get) 840 { 841 uint64_t seed[SZU64]; 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 /* Unscramble the seed. */ 868 scramble_seed (seed, rs->s); 869 870 /* This code now should do correct strides. */ 871 for (size_t i = 0; i < SZ; i++) 872 memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &seed[i], 873 sizeof (GFC_UINTEGER_8)); 874 875 get->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(get, 0)] = rs->p; 876 } 877 878 else 879 { 880 __gthread_mutex_lock (&random_lock); 881 882 /* From the standard: "If no argument is present, the processor assigns 883 a processor-dependent value to the seed." */ 884 if (size == NULL && put == NULL && get == NULL) 885 { 886 master_init = false; 887 init_rand_state (rs, true); 888 } 889 890 if (put != NULL) 891 { 892 /* If the rank of the array is not 1, abort. */ 893 if (GFC_DESCRIPTOR_RANK (put) != 1) 894 runtime_error ("Array rank of PUT is not 1."); 895 896 /* If the array is too small, abort. */ 897 if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ + 1) 898 runtime_error ("Array size of PUT is too small."); 899 900 /* This code now should do correct strides. */ 901 for (size_t i = 0; i < SZ; i++) 902 memcpy (&seed[i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]), 903 sizeof (GFC_UINTEGER_8)); 904 905 scramble_seed (master_state, seed); 906 njumps = 0; 907 master_init = true; 908 init_rand_state (rs, true); 909 rs->p = put->base_addr[SZ * GFC_DESCRIPTOR_STRIDE(put, 0)] & 15; 910 } 911 912 913 __gthread_mutex_unlock (&random_lock); 914 } 915 } 916 iexport(random_seed_i8); 917 918 919 #if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS 920 static void __attribute__((constructor)) 921 constructor_random (void) 922 { 923 #ifndef __GTHREAD_MUTEX_INIT 924 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock); 925 #endif 926 if (__gthread_active_p ()) 927 __gthread_key_create (&rand_state_key, &free); 928 } 929 #endif 930 931 #ifdef __GTHREADS 932 static void __attribute__((destructor)) 933 destructor_random (void) 934 { 935 if (__gthread_active_p ()) 936 __gthread_key_delete (rand_state_key); 937 } 938 #endif 939