xref: /netbsd-src/external/gpl3/gcc/dist/libgfortran/intrinsics/random.c (revision b1e838363e3c6fc78a55519254d99869742dd33c)
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
rnumber_4(GFC_REAL_4 * f,GFC_UINTEGER_4 v)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
rnumber_8(GFC_REAL_8 * f,GFC_UINTEGER_8 v)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
rnumber_10(GFC_REAL_10 * f,GFC_UINTEGER_8 v)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
rnumber_16(GFC_REAL_16 * f,GFC_UINTEGER_8 v1,GFC_UINTEGER_8 v2)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
rnumber_17(GFC_REAL_17 * f,GFC_UINTEGER_8 v1,GFC_UINTEGER_8 v2)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*
get_rand_state(void)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
rotl(const uint64_t x,int k)258 rotl (const uint64_t x, int k)
259 {
260 	return (x << k) | (x >> (64 - k));
261 }
262 
263 
264 static uint64_t
prng_next(prng_state * rs)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
jump(prng_state * rs)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
splitmix64(uint64_t x)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
getosrandom(void * buf,size_t buflen)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
init_rand_state(prng_state * rs,const bool locked)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
random_r4(GFC_REAL_4 * x)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
random_r8(GFC_REAL_8 * x)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
random_r10(GFC_REAL_10 * x)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
random_r16(GFC_REAL_16 * x)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
random_r17(GFC_REAL_17 * x)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
arandom_r4(gfc_array_r4 * x)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
arandom_r8(gfc_array_r8 * x)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
arandom_r10(gfc_array_r10 * x)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
arandom_r16(gfc_array_r16 * x)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
arandom_r17(gfc_array_r17 * x)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
scramble_seed(uint64_t * dest,const uint64_t * src)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
random_seed_i4(GFC_INTEGER_4 * size,gfc_array_i4 * put,gfc_array_i4 * get)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
random_seed_i8(GFC_INTEGER_8 * size,gfc_array_i8 * put,gfc_array_i8 * get)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))
constructor_random(void)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))
destructor_random(void)1048 destructor_random (void)
1049 {
1050   if (__gthread_active_p ())
1051     __gthread_key_delete (rand_state_key);
1052 }
1053 #endif
1054