1 /* Implementation of the RANDOM intrinsics
2 Copyright (C) 2002-2017 Free Software Foundation, Inc.
3 Contributed by Lars Segerlund <seger@linuxmail.org>,
4 Steve Kargl and Janne Blomqvist.
6 This file is part of the GNU Fortran runtime library (libgfortran).
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.
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.
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.
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/>. */
30 #include "libgfortran.h"
44 #include <_mingw.h> /* For __MINGW64_VERSION_MAJOR */
47 extern void random_r4 (GFC_REAL_4
*);
48 iexport_proto(random_r4
);
50 extern void random_r8 (GFC_REAL_8
*);
51 iexport_proto(random_r8
);
53 extern void arandom_r4 (gfc_array_r4
*);
54 export_proto(arandom_r4
);
56 extern void arandom_r8 (gfc_array_r8
*);
57 export_proto(arandom_r8
);
59 #ifdef HAVE_GFC_REAL_10
61 extern void random_r10 (GFC_REAL_10
*);
62 iexport_proto(random_r10
);
64 extern void arandom_r10 (gfc_array_r10
*);
65 export_proto(arandom_r10
);
69 #ifdef HAVE_GFC_REAL_16
71 extern void random_r16 (GFC_REAL_16
*);
72 iexport_proto(random_r16
);
74 extern void arandom_r16 (gfc_array_r16
*);
75 export_proto(arandom_r16
);
79 #ifdef __GTHREAD_MUTEX_INIT
80 static __gthread_mutex_t random_lock
= __GTHREAD_MUTEX_INIT
;
82 static __gthread_mutex_t random_lock
;
85 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
86 GFC_REAL_* types in the range of [0,1). If GFC_REAL_*_RADIX are 2
87 or 16, respectively, we mask off the bits that don't fit into the
88 correct GFC_REAL_*, convert to the real type, then multiply by the
93 rnumber_4 (GFC_REAL_4
*f
, GFC_UINTEGER_4 v
)
96 #if GFC_REAL_4_RADIX == 2
97 mask
= ~ (GFC_UINTEGER_4
) 0u << (32 - GFC_REAL_4_DIGITS
);
98 #elif GFC_REAL_4_RADIX == 16
99 mask
= ~ (GFC_UINTEGER_4
) 0u << ((8 - GFC_REAL_4_DIGITS
) * 4);
101 #error "GFC_REAL_4_RADIX has unknown value"
104 *f
= (GFC_REAL_4
) v
* GFC_REAL_4_LITERAL(0x1.p
-32);
108 rnumber_8 (GFC_REAL_8
*f
, GFC_UINTEGER_8 v
)
111 #if GFC_REAL_8_RADIX == 2
112 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_8_DIGITS
);
113 #elif GFC_REAL_8_RADIX == 16
114 mask
= ~ (GFC_UINTEGER_8
) 0u << (16 - GFC_REAL_8_DIGITS
) * 4);
116 #error "GFC_REAL_8_RADIX has unknown value"
119 *f
= (GFC_REAL_8
) v
* GFC_REAL_8_LITERAL(0x1.p
-64);
122 #ifdef HAVE_GFC_REAL_10
125 rnumber_10 (GFC_REAL_10
*f
, GFC_UINTEGER_8 v
)
128 #if GFC_REAL_10_RADIX == 2
129 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_10_DIGITS
);
130 #elif GFC_REAL_10_RADIX == 16
131 mask
= ~ (GFC_UINTEGER_10
) 0u << ((16 - GFC_REAL_10_DIGITS
) * 4);
133 #error "GFC_REAL_10_RADIX has unknown value"
136 *f
= (GFC_REAL_10
) v
* GFC_REAL_10_LITERAL(0x1.p
-64);
140 #ifdef HAVE_GFC_REAL_16
142 /* For REAL(KIND=16), we only need to mask off the lower bits. */
145 rnumber_16 (GFC_REAL_16
*f
, GFC_UINTEGER_8 v1
, GFC_UINTEGER_8 v2
)
148 #if GFC_REAL_16_RADIX == 2
149 mask
= ~ (GFC_UINTEGER_8
) 0u << (128 - GFC_REAL_16_DIGITS
);
150 #elif GFC_REAL_16_RADIX == 16
151 mask
= ~ (GFC_UINTEGER_8
) 0u << ((32 - GFC_REAL_16_DIGITS
) * 4);
153 #error "GFC_REAL_16_RADIX has unknown value"
156 *f
= (GFC_REAL_16
) v1
* GFC_REAL_16_LITERAL(0x1.p
-64)
157 + (GFC_REAL_16
) v2
* GFC_REAL_16_LITERAL(0x1.p
-128);
164 We use the xorshift1024* generator, a fast high-quality generator
167 - passes TestU1 without any failures
169 - provides a "jump" function making it easy to provide many
170 independent parallel streams.
172 - Long period of 2**1024 - 1
174 A description can be found at
176 http://vigna.di.unimi.it/ftp/papers/xorshift.pdf
180 http://arxiv.org/abs/1402.6246
182 The paper includes public domain source code which is the basis for
183 the implementation below.
192 xorshift1024star_state
;
195 /* master_init, njumps, and master_state are the only variables
196 protected by random_lock. */
197 static bool master_init
;
198 static unsigned njumps
; /* How many times we have jumped. */
199 static uint64_t master_state
[] = {
200 0xad63fa1ed3b55f36ULL
, 0xd94473e78978b497ULL
, 0xbc60592a98172477ULL
,
201 0xa3de7c6e81265301ULL
, 0x586640c5e785af27ULL
, 0x7a2a3f63b67ce5eaULL
,
202 0x9fde969f922d9b82ULL
, 0xe6fe34379b3f3822ULL
, 0x6c277eac3e99b6c2ULL
,
203 0x9197290ab0d3f069ULL
, 0xdb227302f6c25576ULL
, 0xee0209aee527fae9ULL
,
204 0x675666a793cd05b9ULL
, 0xd048c99fbc70c20fULL
, 0x775f8c3dba385ef5ULL
,
205 0x625288bc262faf33ULL
209 static __gthread_key_t rand_state_key
;
211 static xorshift1024star_state
*
212 get_rand_state (void)
214 /* For single threaded apps. */
215 static xorshift1024star_state rand_state
;
217 if (__gthread_active_p ())
219 void* p
= __gthread_getspecific (rand_state_key
);
222 p
= xcalloc (1, sizeof (xorshift1024star_state
));
223 __gthread_setspecific (rand_state_key
, p
);
233 xorshift1024star (xorshift1024star_state
* rs
)
236 const uint64_t s0
= rs
->s
[p
];
237 uint64_t s1
= rs
->s
[p
= (p
+ 1) & 15];
239 rs
->s
[p
] = s1
^ s0
^ (s1
>> 11) ^ (s0
>> 30);
241 return rs
->s
[p
] * UINT64_C(1181783497276652981);
245 /* This is the jump function for the generator. It is equivalent to
246 2^512 calls to xorshift1024star(); it can be used to generate 2^512
247 non-overlapping subsequences for parallel computations. */
250 jump (xorshift1024star_state
* rs
)
252 static const uint64_t JUMP
[] = {
253 0x84242f96eca9c41dULL
, 0xa3c65b8776f96855ULL
, 0x5b34a39f070b5837ULL
,
254 0x4489affce4f31a1eULL
, 0x2ffeeb0a48316f40ULL
, 0xdc2d9891fe68c022ULL
,
255 0x3659132bb12fea70ULL
, 0xaac17d8efa43cab8ULL
, 0xc4cb815590989b13ULL
,
256 0x5ee975283d71c93bULL
, 0x691548c86c1bd540ULL
, 0x7910c41d10a1e6a5ULL
,
257 0x0b5fc64563b3e2a8ULL
, 0x047f7684e9fc949dULL
, 0xb99181f2d8f685caULL
,
258 0x284600e3f30e38c3ULL
261 uint64_t t
[16] = { 0 };
262 for(unsigned int i
= 0; i
< sizeof JUMP
/ sizeof *JUMP
; i
++)
263 for(int b
= 0; b
< 64; b
++)
265 if (JUMP
[i
] & 1ULL << b
)
266 for(int j
= 0; j
< 16; j
++)
267 t
[j
] ^= rs
->s
[(j
+ rs
->p
) & 15];
268 xorshift1024star (rs
);
270 for(int j
= 0; j
< 16; j
++)
271 rs
->s
[(j
+ rs
->p
) & 15] = t
[j
];
275 /* Super-simple LCG generator used in getosrandom () if /dev/urandom
278 #define M 2147483647 /* 2^31 - 1 (A large prime number) */
279 #define A 16807 /* Prime root of M, passes statistical tests and produces a full cycle */
280 #define Q 127773 /* M / A (To avoid overflow on A * seed) */
281 #define R 2836 /* M % A (To avoid overflow on A * seed) */
283 __attribute__((unused
)) static uint32_t
284 lcg_parkmiller(uint32_t seed
)
286 uint32_t hi
= seed
/ Q
;
287 uint32_t lo
= seed
% Q
;
288 int32_t test
= A
* lo
- R
* hi
;
300 /* Get some random bytes from the operating system in order to seed
304 getosrandom (void *buf
, size_t buflen
)
306 /* rand_s is available in MinGW-w64 but not plain MinGW. */
307 #if defined(__MINGW64_VERSION_MAJOR)
308 unsigned int* b
= buf
;
309 for (unsigned i
= 0; i
< buflen
/ sizeof (unsigned int); i
++)
314 TODO: When glibc adds a wrapper for the getrandom() system call
315 on Linux, one could use that.
317 TODO: One could use getentropy() on OpenBSD. */
318 int flags
= O_RDONLY
;
322 int fd
= open("/dev/urandom", flags
);
325 int res
= read(fd
, buf
, buflen
);
329 uint32_t seed
= 1234567890;
332 if (gf_gettime (&secs
, &usecs
) == 0)
338 pid_t pid
= getpid();
342 for (size_t i
= 0; i
< buflen
/ sizeof (uint32_t); i
++)
345 seed
= lcg_parkmiller (seed
);
348 #endif /* __MINGW64_VERSION_MAJOR */
352 /* Initialize the random number generator for the current thread,
353 using the master state and the number of times we must jump. */
356 init_rand_state (xorshift1024star_state
* rs
, const bool locked
)
359 __gthread_mutex_lock (&random_lock
);
362 getosrandom (master_state
, sizeof (master_state
));
366 memcpy (&rs
->s
, master_state
, sizeof (master_state
));
367 unsigned n
= njumps
++;
369 __gthread_mutex_unlock (&random_lock
);
370 for (unsigned i
= 0; i
< n
; i
++)
376 /* This function produces a REAL(4) value from the uniform distribution
380 random_r4 (GFC_REAL_4
*x
)
382 xorshift1024star_state
* rs
= get_rand_state();
384 if (unlikely (!rs
->init
))
385 init_rand_state (rs
, false);
386 uint64_t r
= xorshift1024star (rs
);
387 /* Take the higher bits, ensuring that a stream of real(4), real(8),
388 and real(10) will be identical (except for precision). */
389 uint32_t high
= (uint32_t) (r
>> 32);
394 /* This function produces a REAL(8) value from the uniform distribution
398 random_r8 (GFC_REAL_8
*x
)
401 xorshift1024star_state
* rs
= get_rand_state();
403 if (unlikely (!rs
->init
))
404 init_rand_state (rs
, false);
405 r
= xorshift1024star (rs
);
410 #ifdef HAVE_GFC_REAL_10
412 /* This function produces a REAL(10) value from the uniform distribution
416 random_r10 (GFC_REAL_10
*x
)
419 xorshift1024star_state
* rs
= get_rand_state();
421 if (unlikely (!rs
->init
))
422 init_rand_state (rs
, false);
423 r
= xorshift1024star (rs
);
430 /* This function produces a REAL(16) value from the uniform distribution
433 #ifdef HAVE_GFC_REAL_16
436 random_r16 (GFC_REAL_16
*x
)
438 GFC_UINTEGER_8 r1
, r2
;
439 xorshift1024star_state
* rs
= get_rand_state();
441 if (unlikely (!rs
->init
))
442 init_rand_state (rs
, false);
443 r1
= xorshift1024star (rs
);
444 r2
= xorshift1024star (rs
);
445 rnumber_16 (x
, r1
, r2
);
452 /* This function fills a REAL(4) array with values from the uniform
453 distribution with range [0,1). */
456 arandom_r4 (gfc_array_r4
*x
)
458 index_type count
[GFC_MAX_DIMENSIONS
];
459 index_type extent
[GFC_MAX_DIMENSIONS
];
460 index_type stride
[GFC_MAX_DIMENSIONS
];
464 xorshift1024star_state
* rs
= get_rand_state();
470 dim
= GFC_DESCRIPTOR_RANK (x
);
472 for (n
= 0; n
< dim
; n
++)
475 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
476 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
483 if (unlikely (!rs
->init
))
484 init_rand_state (rs
, false);
488 /* random_r4 (dest); */
489 uint64_t r
= xorshift1024star (rs
);
490 uint32_t high
= (uint32_t) (r
>> 32);
491 rnumber_4 (dest
, high
);
493 /* Advance to the next element. */
496 /* Advance to the next source element. */
498 while (count
[n
] == extent
[n
])
500 /* When we get to the end of a dimension, reset it and increment
501 the next dimension. */
503 /* We could precalculate these products, but this is a less
504 frequently used path so probably not worth it. */
505 dest
-= stride
[n
] * extent
[n
];
521 /* This function fills a REAL(8) array with values from the uniform
522 distribution with range [0,1). */
525 arandom_r8 (gfc_array_r8
*x
)
527 index_type count
[GFC_MAX_DIMENSIONS
];
528 index_type extent
[GFC_MAX_DIMENSIONS
];
529 index_type stride
[GFC_MAX_DIMENSIONS
];
533 xorshift1024star_state
* rs
= get_rand_state();
538 dim
= GFC_DESCRIPTOR_RANK (x
);
540 for (n
= 0; n
< dim
; n
++)
543 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
544 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
551 if (unlikely (!rs
->init
))
552 init_rand_state (rs
, false);
556 /* random_r8 (dest); */
557 uint64_t r
= xorshift1024star (rs
);
560 /* Advance to the next element. */
563 /* Advance to the next source element. */
565 while (count
[n
] == extent
[n
])
567 /* When we get to the end of a dimension, reset it and increment
568 the next dimension. */
570 /* We could precalculate these products, but this is a less
571 frequently used path so probably not worth it. */
572 dest
-= stride
[n
] * extent
[n
];
588 #ifdef HAVE_GFC_REAL_10
590 /* This function fills a REAL(10) array with values from the uniform
591 distribution with range [0,1). */
594 arandom_r10 (gfc_array_r10
*x
)
596 index_type count
[GFC_MAX_DIMENSIONS
];
597 index_type extent
[GFC_MAX_DIMENSIONS
];
598 index_type stride
[GFC_MAX_DIMENSIONS
];
602 xorshift1024star_state
* rs
= get_rand_state();
607 dim
= GFC_DESCRIPTOR_RANK (x
);
609 for (n
= 0; n
< dim
; n
++)
612 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
613 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
620 if (unlikely (!rs
->init
))
621 init_rand_state (rs
, false);
625 /* random_r10 (dest); */
626 uint64_t r
= xorshift1024star (rs
);
627 rnumber_10 (dest
, r
);
629 /* Advance to the next element. */
632 /* Advance to the next source element. */
634 while (count
[n
] == extent
[n
])
636 /* When we get to the end of a dimension, reset it and increment
637 the next dimension. */
639 /* We could precalculate these products, but this is a less
640 frequently used path so probably not worth it. */
641 dest
-= stride
[n
] * extent
[n
];
659 #ifdef HAVE_GFC_REAL_16
661 /* This function fills a REAL(16) array with values from the uniform
662 distribution with range [0,1). */
665 arandom_r16 (gfc_array_r16
*x
)
667 index_type count
[GFC_MAX_DIMENSIONS
];
668 index_type extent
[GFC_MAX_DIMENSIONS
];
669 index_type stride
[GFC_MAX_DIMENSIONS
];
673 xorshift1024star_state
* rs
= get_rand_state();
678 dim
= GFC_DESCRIPTOR_RANK (x
);
680 for (n
= 0; n
< dim
; n
++)
683 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
684 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
691 if (unlikely (!rs
->init
))
692 init_rand_state (rs
, false);
696 /* random_r16 (dest); */
697 uint64_t r1
= xorshift1024star (rs
);
698 uint64_t r2
= xorshift1024star (rs
);
699 rnumber_16 (dest
, r1
, r2
);
701 /* Advance to the next element. */
704 /* Advance to the next source element. */
706 while (count
[n
] == extent
[n
])
708 /* When we get to the end of a dimension, reset it and increment
709 the next dimension. */
711 /* We could precalculate these products, but this is a less
712 frequently used path so probably not worth it. */
713 dest
-= stride
[n
] * extent
[n
];
732 /* Number of elements in master_state array. */
733 #define SZU64 (sizeof (master_state) / sizeof (uint64_t))
736 /* Keys for scrambling the seed in order to avoid poor seeds. */
738 static const uint64_t xor_keys
[] = {
739 0xbd0c5b6e50c2df49ULL
, 0xd46061cd46e1df38ULL
, 0xbb4f4d4ed6103544ULL
,
740 0x114a583d0756ad39ULL
, 0x4b5ad8623d0aaab6ULL
, 0x3f2ed7afbe0c0f21ULL
,
741 0xdec83fd65f113445ULL
, 0x3824f8fbc4f10d24ULL
, 0x5d9025af05878911ULL
,
742 0x500bc46b540340e9ULL
, 0x8bd53298e0d00530ULL
, 0x57886e40a952e06aULL
,
743 0x926e76c88e31cdb6ULL
, 0xbd0724dac0a3a5f9ULL
, 0xc5c8981b858ab796ULL
,
744 0xbb12ab2694c2b32cULL
748 /* Since a XOR cipher is symmetric, we need only one routine, and we
749 can use it both for encryption and decryption. */
752 scramble_seed (uint64_t *dest
, const uint64_t *src
)
754 for (int i
= 0; i
< (int) SZU64
; i
++)
755 dest
[i
] = src
[i
] ^ xor_keys
[i
];
759 /* random_seed is used to seed the PRNG with either a default
760 set of seeds or user specified set of seeds. random_seed
761 must be called with no argument or exactly one argument. */
764 random_seed_i4 (GFC_INTEGER_4
*size
, gfc_array_i4
*put
, gfc_array_i4
*get
)
766 uint64_t seed
[SZU64
];
767 #define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_4))
769 /* Check that we only have one argument present. */
770 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
771 runtime_error ("RANDOM_SEED should have at most one argument present.");
776 xorshift1024star_state
* rs
= get_rand_state();
778 /* Return the seed to GET data. */
781 /* If the rank of the array is not 1, abort. */
782 if (GFC_DESCRIPTOR_RANK (get
) != 1)
783 runtime_error ("Array rank of GET is not 1.");
785 /* If the array is too small, abort. */
786 if (GFC_DESCRIPTOR_EXTENT(get
,0) < (index_type
) SZ
+ 1)
787 runtime_error ("Array size of GET is too small.");
790 init_rand_state (rs
, false);
792 /* Unscramble the seed. */
793 scramble_seed (seed
, rs
->s
);
795 /* Then copy it back to the user variable. */
796 for (size_t i
= 0; i
< SZ
; i
++)
797 memcpy (&(get
->base_addr
[(SZ
- 1 - i
) * GFC_DESCRIPTOR_STRIDE(get
,0)]),
798 (unsigned char*) seed
+ i
* sizeof(GFC_UINTEGER_4
),
799 sizeof(GFC_UINTEGER_4
));
801 /* Finally copy the value of p after the seed. */
802 get
->base_addr
[SZ
* GFC_DESCRIPTOR_STRIDE(get
, 0)] = rs
->p
;
807 __gthread_mutex_lock (&random_lock
);
809 /* From the standard: "If no argument is present, the processor assigns
810 a processor-dependent value to the seed." */
811 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
814 init_rand_state (rs
, true);
819 /* If the rank of the array is not 1, abort. */
820 if (GFC_DESCRIPTOR_RANK (put
) != 1)
821 runtime_error ("Array rank of PUT is not 1.");
823 /* If the array is too small, abort. */
824 if (GFC_DESCRIPTOR_EXTENT(put
,0) < (index_type
) SZ
+ 1)
825 runtime_error ("Array size of PUT is too small.");
827 /* We copy the seed given by the user. */
828 for (size_t i
= 0; i
< SZ
; i
++)
829 memcpy ((unsigned char*) seed
+ i
* sizeof(GFC_UINTEGER_4
),
830 &(put
->base_addr
[(SZ
- 1 - i
) * GFC_DESCRIPTOR_STRIDE(put
,0)]),
831 sizeof(GFC_UINTEGER_4
));
833 /* We put it after scrambling the bytes, to paper around users who
834 provide seeds with quality only in the lower or upper part. */
835 scramble_seed (master_state
, seed
);
838 init_rand_state (rs
, true);
840 rs
->p
= put
->base_addr
[SZ
* GFC_DESCRIPTOR_STRIDE(put
, 0)] & 15;
843 __gthread_mutex_unlock (&random_lock
);
847 iexport(random_seed_i4
);
851 random_seed_i8 (GFC_INTEGER_8
*size
, gfc_array_i8
*put
, gfc_array_i8
*get
)
853 uint64_t seed
[SZU64
];
855 /* Check that we only have one argument present. */
856 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
857 runtime_error ("RANDOM_SEED should have at most one argument present.");
859 #define SZ (sizeof (master_state) / sizeof (GFC_INTEGER_8))
863 xorshift1024star_state
* rs
= get_rand_state();
865 /* Return the seed to GET data. */
868 /* If the rank of the array is not 1, abort. */
869 if (GFC_DESCRIPTOR_RANK (get
) != 1)
870 runtime_error ("Array rank of GET is not 1.");
872 /* If the array is too small, abort. */
873 if (GFC_DESCRIPTOR_EXTENT(get
,0) < (index_type
) SZ
+ 1)
874 runtime_error ("Array size of GET is too small.");
877 init_rand_state (rs
, false);
879 /* Unscramble the seed. */
880 scramble_seed (seed
, rs
->s
);
882 /* This code now should do correct strides. */
883 for (size_t i
= 0; i
< SZ
; i
++)
884 memcpy (&(get
->base_addr
[i
* GFC_DESCRIPTOR_STRIDE(get
,0)]), &seed
[i
],
885 sizeof (GFC_UINTEGER_8
));
887 get
->base_addr
[SZ
* GFC_DESCRIPTOR_STRIDE(get
, 0)] = rs
->p
;
892 __gthread_mutex_lock (&random_lock
);
894 /* From the standard: "If no argument is present, the processor assigns
895 a processor-dependent value to the seed." */
896 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
899 init_rand_state (rs
, true);
904 /* If the rank of the array is not 1, abort. */
905 if (GFC_DESCRIPTOR_RANK (put
) != 1)
906 runtime_error ("Array rank of PUT is not 1.");
908 /* If the array is too small, abort. */
909 if (GFC_DESCRIPTOR_EXTENT(put
,0) < (index_type
) SZ
+ 1)
910 runtime_error ("Array size of PUT is too small.");
912 /* This code now should do correct strides. */
913 for (size_t i
= 0; i
< SZ
; i
++)
914 memcpy (&seed
[i
], &(put
->base_addr
[i
* GFC_DESCRIPTOR_STRIDE(put
,0)]),
915 sizeof (GFC_UINTEGER_8
));
917 scramble_seed (master_state
, seed
);
920 init_rand_state (rs
, true);
921 rs
->p
= put
->base_addr
[SZ
* GFC_DESCRIPTOR_STRIDE(put
, 0)] & 15;
925 __gthread_mutex_unlock (&random_lock
);
928 iexport(random_seed_i8
);
931 #if !defined __GTHREAD_MUTEX_INIT || defined __GTHREADS
932 static void __attribute__((constructor
))
933 constructor_random (void)
935 #ifndef __GTHREAD_MUTEX_INIT
936 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock
);
938 if (__gthread_active_p ())
939 __gthread_key_create (&rand_state_key
, &free
);
944 static void __attribute__((destructor
))
945 destructor_random (void)
947 if (__gthread_active_p ())
948 __gthread_key_delete (rand_state_key
);