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
rnumber_4(GFC_REAL_4 * f,GFC_UINTEGER_4 v)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
rnumber_8(GFC_REAL_8 * f,GFC_UINTEGER_8 v)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
rnumber_10(GFC_REAL_10 * f,GFC_UINTEGER_8 v)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
rnumber_16(GFC_REAL_16 * f,GFC_UINTEGER_8 v1,GFC_UINTEGER_8 v2)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*
get_rand_state(void)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
rotl(const uint64_t x,int k)227 rotl (const uint64_t x, int k)
228 {
229 return (x << k) | (x >> (64 - k));
230 }
231
232
233 static uint64_t
prng_next(prng_state * rs)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
jump(prng_state * rs)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
splitmix64(uint64_t x)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
getosrandom(void * buf,size_t buflen)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
init_rand_state(prng_state * rs,const bool locked)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
random_r4(GFC_REAL_4 * x)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
random_r8(GFC_REAL_8 * x)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
random_r10(GFC_REAL_10 * x)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
random_r16(GFC_REAL_16 * x)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
arandom_r4(gfc_array_r4 * x)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
arandom_r8(gfc_array_r8 * x)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
arandom_r10(gfc_array_r10 * x)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
arandom_r16(gfc_array_r16 * x)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
scramble_seed(uint64_t * dest,const uint64_t * src)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
random_seed_i4(GFC_INTEGER_4 * size,gfc_array_i4 * put,gfc_array_i4 * get)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
random_seed_i8(GFC_INTEGER_8 * size,gfc_array_i8 * put,gfc_array_i8 * get)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))
constructor_random(void)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))
destructor_random(void)922 destructor_random (void)
923 {
924 if (__gthread_active_p ())
925 __gthread_key_delete (rand_state_key);
926 }
927 #endif
928