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