1 /* Implementation of the RANDOM intrinsics 2 Copyright (C) 2002-2020 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 727 /* Keys for scrambling the seed in order to avoid poor seeds. */ 728 729 static const uint64_t xor_keys[] = { 730 0xbd0c5b6e50c2df49ULL, 0xd46061cd46e1df38ULL, 0xbb4f4d4ed6103544ULL, 731 0x114a583d0756ad39ULL 732 }; 733 734 735 /* Since a XOR cipher is symmetric, we need only one routine, and we 736 can use it both for encryption and decryption. */ 737 738 static void 739 scramble_seed (uint64_t *dest, const uint64_t *src) 740 { 741 for (size_t i = 0; i < SZU64; i++) 742 dest[i] = src[i] ^ xor_keys[i]; 743 } 744 745 746 /* random_seed is used to seed the PRNG with either a default 747 set of seeds or user specified set of seeds. random_seed 748 must be called with no argument or exactly one argument. */ 749 750 void 751 random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) 752 { 753 uint64_t seed[SZU64]; 754 #define SZ (sizeof (master_state.s) / sizeof (GFC_INTEGER_4)) 755 756 /* Check that we only have one argument present. */ 757 if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) 758 runtime_error ("RANDOM_SEED should have at most one argument present."); 759 760 if (size != NULL) 761 *size = SZ; 762 763 prng_state* rs = get_rand_state(); 764 765 /* Return the seed to GET data. */ 766 if (get != NULL) 767 { 768 /* If the rank of the array is not 1, abort. */ 769 if (GFC_DESCRIPTOR_RANK (get) != 1) 770 runtime_error ("Array rank of GET is not 1."); 771 772 /* If the array is too small, abort. */ 773 if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ) 774 runtime_error ("Array size of GET is too small."); 775 776 if (!rs->init) 777 init_rand_state (rs, false); 778 779 /* Unscramble the seed. */ 780 scramble_seed (seed, rs->s); 781 782 /* Then copy it back to the user variable. */ 783 for (size_t i = 0; i < SZ ; i++) 784 memcpy (&(get->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]), 785 (unsigned char*) seed + i * sizeof(GFC_UINTEGER_4), 786 sizeof(GFC_UINTEGER_4)); 787 } 788 789 else 790 { 791 __gthread_mutex_lock (&random_lock); 792 793 /* From the standard: "If no argument is present, the processor assigns 794 a processor-dependent value to the seed." */ 795 if (size == NULL && put == NULL && get == NULL) 796 { 797 master_state.init = false; 798 init_rand_state (rs, true); 799 } 800 801 if (put != NULL) 802 { 803 /* If the rank of the array is not 1, abort. */ 804 if (GFC_DESCRIPTOR_RANK (put) != 1) 805 runtime_error ("Array rank of PUT is not 1."); 806 807 /* If the array is too small, abort. */ 808 if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ) 809 runtime_error ("Array size of PUT is too small."); 810 811 /* We copy the seed given by the user. */ 812 for (size_t i = 0; i < SZ; i++) 813 memcpy ((unsigned char*) seed + i * sizeof(GFC_UINTEGER_4), 814 &(put->base_addr[(SZ - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]), 815 sizeof(GFC_UINTEGER_4)); 816 817 /* We put it after scrambling the bytes, to paper around users who 818 provide seeds with quality only in the lower or upper part. */ 819 scramble_seed (master_state.s, seed); 820 master_state.init = true; 821 init_rand_state (rs, true); 822 } 823 824 __gthread_mutex_unlock (&random_lock); 825 } 826 #undef SZ 827 } 828 iexport(random_seed_i4); 829 830 831 void 832 random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get) 833 { 834 uint64_t seed[SZU64]; 835 836 /* Check that we only have one argument present. */ 837 if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) 838 runtime_error ("RANDOM_SEED should have at most one argument present."); 839 840 #define SZ (sizeof (master_state.s) / sizeof (GFC_INTEGER_8)) 841 if (size != NULL) 842 *size = SZ; 843 844 prng_state* rs = get_rand_state(); 845 846 /* Return the seed to GET data. */ 847 if (get != NULL) 848 { 849 /* If the rank of the array is not 1, abort. */ 850 if (GFC_DESCRIPTOR_RANK (get) != 1) 851 runtime_error ("Array rank of GET is not 1."); 852 853 /* If the array is too small, abort. */ 854 if (GFC_DESCRIPTOR_EXTENT(get,0) < (index_type) SZ) 855 runtime_error ("Array size of GET is too small."); 856 857 if (!rs->init) 858 init_rand_state (rs, false); 859 860 /* Unscramble the seed. */ 861 scramble_seed (seed, rs->s); 862 863 /* This code now should do correct strides. */ 864 for (size_t i = 0; i < SZ; i++) 865 memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &seed[i], 866 sizeof (GFC_UINTEGER_8)); 867 } 868 869 else 870 { 871 __gthread_mutex_lock (&random_lock); 872 873 /* From the standard: "If no argument is present, the processor assigns 874 a processor-dependent value to the seed." */ 875 if (size == NULL && put == NULL && get == NULL) 876 { 877 master_state.init = false; 878 init_rand_state (rs, true); 879 } 880 881 if (put != NULL) 882 { 883 /* If the rank of the array is not 1, abort. */ 884 if (GFC_DESCRIPTOR_RANK (put) != 1) 885 runtime_error ("Array rank of PUT is not 1."); 886 887 /* If the array is too small, abort. */ 888 if (GFC_DESCRIPTOR_EXTENT(put,0) < (index_type) SZ) 889 runtime_error ("Array size of PUT is too small."); 890 891 /* This code now should do correct strides. */ 892 for (size_t i = 0; i < SZ; i++) 893 memcpy (&seed[i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]), 894 sizeof (GFC_UINTEGER_8)); 895 896 scramble_seed (master_state.s, seed); 897 master_state.init = true; 898 init_rand_state (rs, true); 899 } 900 901 902 __gthread_mutex_unlock (&random_lock); 903 } 904 } 905 iexport(random_seed_i8); 906 907 908 #if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS 909 static void __attribute__((constructor)) 910 constructor_random (void) 911 { 912 #ifndef __GTHREAD_MUTEX_INIT 913 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock); 914 #endif 915 if (__gthread_active_p ()) 916 __gthread_key_create (&rand_state_key, &free); 917 } 918 #endif 919 920 #ifdef __GTHREADS 921 static void __attribute__((destructor)) 922 destructor_random (void) 923 { 924 if (__gthread_active_p ()) 925 __gthread_key_delete (rand_state_key); 926 } 927 #endif 928