1 /* Implementation of the RANDOM intrinsics
2 Copyright 2002, 2004, 2005, 2006, 2007, 2009, 2010
3 Free Software Foundation, Inc.
4 Contributed by Lars Segerlund <seger@linuxmail.org>
7 This file is part of the GNU Fortran 95 runtime library (libgfortran).
9 Libgfortran is free software; you can redistribute it and/or
10 modify it under the terms of the GNU General Public
11 License as published by the Free Software Foundation; either
12 version 3 of the License, or (at your option) any later version.
14 Ligbfortran is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
19 Under Section 7 of GPL version 3, you are granted additional
20 permissions described in the GCC Runtime Library Exception, version
21 3.1, as published by the Free Software Foundation.
23 You should have received a copy of the GNU General Public License and
24 a copy of the GCC Runtime Library Exception along with this program;
25 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
26 <http://www.gnu.org/licenses/>. */
28 #include "libgfortran.h"
32 extern void random_r4 (GFC_REAL_4
*);
33 iexport_proto(random_r4
);
35 extern void random_r8 (GFC_REAL_8
*);
36 iexport_proto(random_r8
);
38 extern void arandom_r4 (gfc_array_r4
*);
39 export_proto(arandom_r4
);
41 extern void arandom_r8 (gfc_array_r8
*);
42 export_proto(arandom_r8
);
44 #ifdef HAVE_GFC_REAL_10
46 extern void random_r10 (GFC_REAL_10
*);
47 iexport_proto(random_r10
);
49 extern void arandom_r10 (gfc_array_r10
*);
50 export_proto(arandom_r10
);
54 #ifdef HAVE_GFC_REAL_16
56 extern void random_r16 (GFC_REAL_16
*);
57 iexport_proto(random_r16
);
59 extern void arandom_r16 (gfc_array_r16
*);
60 export_proto(arandom_r16
);
64 #ifdef __GTHREAD_MUTEX_INIT
65 static __gthread_mutex_t random_lock
= __GTHREAD_MUTEX_INIT
;
67 static __gthread_mutex_t random_lock
;
70 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
71 GFC_REAL_* types in the range of [0,1). If GFC_REAL_*_RADIX are 2
72 or 16, respectively, we mask off the bits that don't fit into the
73 correct GFC_REAL_*, convert to the real type, then multiply by the
78 rnumber_4 (GFC_REAL_4
*f
, GFC_UINTEGER_4 v
)
81 #if GFC_REAL_4_RADIX == 2
82 mask
= ~ (GFC_UINTEGER_4
) 0u << (32 - GFC_REAL_4_DIGITS
);
83 #elif GFC_REAL_4_RADIX == 16
84 mask
= ~ (GFC_UINTEGER_4
) 0u << ((8 - GFC_REAL_4_DIGITS
) * 4);
86 #error "GFC_REAL_4_RADIX has unknown value"
89 *f
= (GFC_REAL_4
) v
* GFC_REAL_4_LITERAL(0x1.p
-32);
93 rnumber_8 (GFC_REAL_8
*f
, GFC_UINTEGER_8 v
)
96 #if GFC_REAL_8_RADIX == 2
97 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_8_DIGITS
);
98 #elif GFC_REAL_8_RADIX == 16
99 mask
= ~ (GFC_UINTEGER_8
) 0u << (16 - GFC_REAL_8_DIGITS
) * 4);
101 #error "GFC_REAL_8_RADIX has unknown value"
104 *f
= (GFC_REAL_8
) v
* GFC_REAL_8_LITERAL(0x1.p
-64);
107 #ifdef HAVE_GFC_REAL_10
110 rnumber_10 (GFC_REAL_10
*f
, GFC_UINTEGER_8 v
)
113 #if GFC_REAL_10_RADIX == 2
114 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_10_DIGITS
);
115 #elif GFC_REAL_10_RADIX == 16
116 mask
= ~ (GFC_UINTEGER_10
) 0u << ((16 - GFC_REAL_10_DIGITS
) * 4);
118 #error "GFC_REAL_10_RADIX has unknown value"
121 *f
= (GFC_REAL_10
) v
* GFC_REAL_10_LITERAL(0x1.p
-64);
125 #ifdef HAVE_GFC_REAL_16
127 /* For REAL(KIND=16), we only need to mask off the lower bits. */
130 rnumber_16 (GFC_REAL_16
*f
, GFC_UINTEGER_8 v1
, GFC_UINTEGER_8 v2
)
133 #if GFC_REAL_16_RADIX == 2
134 mask
= ~ (GFC_UINTEGER_8
) 0u << (128 - GFC_REAL_16_DIGITS
);
135 #elif GFC_REAL_16_RADIX == 16
136 mask
= ~ (GFC_UINTEGER_8
) 0u << ((32 - GFC_REAL_16_DIGITS
) * 4);
138 #error "GFC_REAL_16_RADIX has unknown value"
141 *f
= (GFC_REAL_16
) v1
* GFC_REAL_16_LITERAL(0x1.p
-64)
142 + (GFC_REAL_16
) v2
* GFC_REAL_16_LITERAL(0x1.p
-128);
145 /* libgfortran previously had a Mersenne Twister, taken from the paper:
147 Mersenne Twister: 623-dimensionally equidistributed
148 uniform pseudorandom generator.
150 by Makoto Matsumoto & Takuji Nishimura
151 which appeared in the: ACM Transactions on Modelling and Computer
152 Simulations: Special Issue on Uniform Random Number
153 Generation. ( Early in 1998 ).
155 The Mersenne Twister code was replaced due to
157 (1) Simple user specified seeds lead to really bad sequences for
158 nearly 100000 random numbers.
159 (2) open(), read(), and close() were not properly declared via header
161 (3) The global index i was abused and caused unexpected behavior with
166 libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
167 random number generator. This PRNG combines:
169 (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
171 (2) A 3-shift shift-register generator with a period of 2^32-1,
172 (3) Two 16-bit multiply-with-carry generators with a period of
173 597273182964842497 > 2^59.
175 The overall period exceeds 2^123.
177 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
179 The above web site has an archive of a newsgroup posting from George
180 Marsaglia with the statement:
182 Subject: Random numbers for C: Improvements.
183 Date: Fri, 15 Jan 1999 11:41:47 -0500
184 From: George Marsaglia <geo@stat.fsu.edu>
185 Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
186 References: <369B5E30.65A55FD1@stat.fsu.edu>
187 Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
190 As I hoped, several suggestions have led to
191 improvements in the code for RNG's I proposed for
192 use in C. (See the thread "Random numbers for C: Some
193 suggestions" in previous postings.) The improved code
196 A question of copyright has also been raised. Unlike
197 DIEHARD, there is no copyright on the code below. You
198 are free to use it in any way you want, but you may
199 wish to acknowledge the source, as a courtesy.
201 "There is no copyright on the code below." included the original
204 /* We use three KISS random number generators, with different
206 As a matter of Quality of Implementation, the random numbers
207 we generate for different REAL kinds, starting from the same
208 seed, are always the same up to the precision of these types.
209 We do this by using three generators with different seeds, the
210 first one always for the most significant bits, the second one
211 for bits 33..64 (if present in the REAL kind), and the third one
212 (called twice) for REAL(16). */
214 #define GFC_SL(k, n) ((k)^((k)<<(n)))
215 #define GFC_SR(k, n) ((k)^((k)>>(n)))
217 /* Reference for the seed:
218 From: "George Marsaglia" <g...@stat.fsu.edu>
220 Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
222 The KISS RNG uses four seeds, x, y, z, c,
223 with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
224 except that the two pairs
225 z=0,c=0 and z=2^32-1,c=698769068
226 should be avoided. */
228 /* Any modifications to the seeds that change kiss_size below need to be
229 reflected in check.c (gfc_check_random_seed) to enable correct
230 compile-time checking of PUT size for the RANDOM_SEED intrinsic. */
232 #define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
233 #define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
234 #ifdef HAVE_GFC_REAL_16
235 #define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
238 static GFC_UINTEGER_4 kiss_seed
[] = {
241 #ifdef HAVE_GFC_REAL_16
246 static GFC_UINTEGER_4 kiss_default_seed
[] = {
249 #ifdef HAVE_GFC_REAL_16
254 static const GFC_INTEGER_4 kiss_size
= sizeof(kiss_seed
)/sizeof(kiss_seed
[0]);
256 static GFC_UINTEGER_4
* const kiss_seed_1
= kiss_seed
;
257 static GFC_UINTEGER_4
* const kiss_seed_2
= kiss_seed
+ 4;
259 #ifdef HAVE_GFC_REAL_16
260 static GFC_UINTEGER_4
* const kiss_seed_3
= kiss_seed
+ 8;
263 /* kiss_random_kernel() returns an integer value in the range of
264 (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
265 should be uniform. */
267 static GFC_UINTEGER_4
268 kiss_random_kernel(GFC_UINTEGER_4
* seed
)
272 seed
[0] = 69069 * seed
[0] + 1327217885;
273 seed
[1] = GFC_SL(GFC_SR(GFC_SL(seed
[1],13),17),5);
274 seed
[2] = 18000 * (seed
[2] & 65535) + (seed
[2] >> 16);
275 seed
[3] = 30903 * (seed
[3] & 65535) + (seed
[3] >> 16);
276 kiss
= seed
[0] + seed
[1] + (seed
[2] << 16) + seed
[3];
281 /* This function produces a REAL(4) value from the uniform distribution
285 random_r4 (GFC_REAL_4
*x
)
289 __gthread_mutex_lock (&random_lock
);
290 kiss
= kiss_random_kernel (kiss_seed_1
);
292 __gthread_mutex_unlock (&random_lock
);
296 /* This function produces a REAL(8) value from the uniform distribution
300 random_r8 (GFC_REAL_8
*x
)
304 __gthread_mutex_lock (&random_lock
);
305 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
306 kiss
+= kiss_random_kernel (kiss_seed_2
);
308 __gthread_mutex_unlock (&random_lock
);
312 #ifdef HAVE_GFC_REAL_10
314 /* This function produces a REAL(10) value from the uniform distribution
318 random_r10 (GFC_REAL_10
*x
)
322 __gthread_mutex_lock (&random_lock
);
323 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
324 kiss
+= kiss_random_kernel (kiss_seed_2
);
325 rnumber_10 (x
, kiss
);
326 __gthread_mutex_unlock (&random_lock
);
332 /* This function produces a REAL(16) value from the uniform distribution
335 #ifdef HAVE_GFC_REAL_16
338 random_r16 (GFC_REAL_16
*x
)
340 GFC_UINTEGER_8 kiss1
, kiss2
;
342 __gthread_mutex_lock (&random_lock
);
343 kiss1
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
344 kiss1
+= kiss_random_kernel (kiss_seed_2
);
346 kiss2
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_3
)) << 32;
347 kiss2
+= kiss_random_kernel (kiss_seed_3
);
349 rnumber_16 (x
, kiss1
, kiss2
);
350 __gthread_mutex_unlock (&random_lock
);
356 /* This function fills a REAL(4) array with values from the uniform
357 distribution with range [0,1). */
360 arandom_r4 (gfc_array_r4
*x
)
362 index_type count
[GFC_MAX_DIMENSIONS
];
363 index_type extent
[GFC_MAX_DIMENSIONS
];
364 index_type stride
[GFC_MAX_DIMENSIONS
];
373 dim
= GFC_DESCRIPTOR_RANK (x
);
375 for (n
= 0; n
< dim
; n
++)
378 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
379 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
386 __gthread_mutex_lock (&random_lock
);
390 /* random_r4 (dest); */
391 kiss
= kiss_random_kernel (kiss_seed_1
);
392 rnumber_4 (dest
, kiss
);
394 /* Advance to the next element. */
397 /* Advance to the next source element. */
399 while (count
[n
] == extent
[n
])
401 /* When we get to the end of a dimension, reset it and increment
402 the next dimension. */
404 /* We could precalculate these products, but this is a less
405 frequently used path so probably not worth it. */
406 dest
-= stride
[n
] * extent
[n
];
420 __gthread_mutex_unlock (&random_lock
);
423 /* This function fills a REAL(8) array with values from the uniform
424 distribution with range [0,1). */
427 arandom_r8 (gfc_array_r8
*x
)
429 index_type count
[GFC_MAX_DIMENSIONS
];
430 index_type extent
[GFC_MAX_DIMENSIONS
];
431 index_type stride
[GFC_MAX_DIMENSIONS
];
440 dim
= GFC_DESCRIPTOR_RANK (x
);
442 for (n
= 0; n
< dim
; n
++)
445 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
446 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
453 __gthread_mutex_lock (&random_lock
);
457 /* random_r8 (dest); */
458 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
459 kiss
+= kiss_random_kernel (kiss_seed_2
);
460 rnumber_8 (dest
, kiss
);
462 /* Advance to the next element. */
465 /* Advance to the next source element. */
467 while (count
[n
] == extent
[n
])
469 /* When we get to the end of a dimension, reset it and increment
470 the next dimension. */
472 /* We could precalculate these products, but this is a less
473 frequently used path so probably not worth it. */
474 dest
-= stride
[n
] * extent
[n
];
488 __gthread_mutex_unlock (&random_lock
);
491 #ifdef HAVE_GFC_REAL_10
493 /* This function fills a REAL(10) array with values from the uniform
494 distribution with range [0,1). */
497 arandom_r10 (gfc_array_r10
*x
)
499 index_type count
[GFC_MAX_DIMENSIONS
];
500 index_type extent
[GFC_MAX_DIMENSIONS
];
501 index_type stride
[GFC_MAX_DIMENSIONS
];
510 dim
= GFC_DESCRIPTOR_RANK (x
);
512 for (n
= 0; n
< dim
; n
++)
515 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
516 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
523 __gthread_mutex_lock (&random_lock
);
527 /* random_r10 (dest); */
528 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
529 kiss
+= kiss_random_kernel (kiss_seed_2
);
530 rnumber_10 (dest
, kiss
);
532 /* Advance to the next element. */
535 /* Advance to the next source element. */
537 while (count
[n
] == extent
[n
])
539 /* When we get to the end of a dimension, reset it and increment
540 the next dimension. */
542 /* We could precalculate these products, but this is a less
543 frequently used path so probably not worth it. */
544 dest
-= stride
[n
] * extent
[n
];
558 __gthread_mutex_unlock (&random_lock
);
563 #ifdef HAVE_GFC_REAL_16
565 /* This function fills a REAL(16) array with values from the uniform
566 distribution with range [0,1). */
569 arandom_r16 (gfc_array_r16
*x
)
571 index_type count
[GFC_MAX_DIMENSIONS
];
572 index_type extent
[GFC_MAX_DIMENSIONS
];
573 index_type stride
[GFC_MAX_DIMENSIONS
];
577 GFC_UINTEGER_8 kiss1
, kiss2
;
582 dim
= GFC_DESCRIPTOR_RANK (x
);
584 for (n
= 0; n
< dim
; n
++)
587 stride
[n
] = GFC_DESCRIPTOR_STRIDE(x
,n
);
588 extent
[n
] = GFC_DESCRIPTOR_EXTENT(x
,n
);
595 __gthread_mutex_lock (&random_lock
);
599 /* random_r16 (dest); */
600 kiss1
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
601 kiss1
+= kiss_random_kernel (kiss_seed_2
);
603 kiss2
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_3
)) << 32;
604 kiss2
+= kiss_random_kernel (kiss_seed_3
);
606 rnumber_16 (dest
, kiss1
, kiss2
);
608 /* Advance to the next element. */
611 /* Advance to the next source element. */
613 while (count
[n
] == extent
[n
])
615 /* When we get to the end of a dimension, reset it and increment
616 the next dimension. */
618 /* We could precalculate these products, but this is a less
619 frequently used path so probably not worth it. */
620 dest
-= stride
[n
] * extent
[n
];
634 __gthread_mutex_unlock (&random_lock
);
642 scramble_seed (unsigned char *dest
, unsigned char *src
, int size
)
646 for (i
= 0; i
< size
; i
++)
647 dest
[(i
% 2) * (size
/ 2) + i
/ 2] = src
[i
];
652 unscramble_seed (unsigned char *dest
, unsigned char *src
, int size
)
656 for (i
= 0; i
< size
; i
++)
657 dest
[i
] = src
[(i
% 2) * (size
/ 2) + i
/ 2];
662 /* random_seed is used to seed the PRNG with either a default
663 set of seeds or user specified set of seeds. random_seed
664 must be called with no argument or exactly one argument. */
667 random_seed_i4 (GFC_INTEGER_4
*size
, gfc_array_i4
*put
, gfc_array_i4
*get
)
670 unsigned char seed
[4*kiss_size
];
672 __gthread_mutex_lock (&random_lock
);
674 /* Check that we only have one argument present. */
675 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
676 runtime_error ("RANDOM_SEED should have at most one argument present.");
678 /* From the standard: "If no argument is present, the processor assigns
679 a processor-dependent value to the seed." */
680 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
681 for (i
= 0; i
< kiss_size
; i
++)
682 kiss_seed
[i
] = kiss_default_seed
[i
];
689 /* If the rank of the array is not 1, abort. */
690 if (GFC_DESCRIPTOR_RANK (put
) != 1)
691 runtime_error ("Array rank of PUT is not 1.");
693 /* If the array is too small, abort. */
694 if (GFC_DESCRIPTOR_EXTENT(put
,0) < kiss_size
)
695 runtime_error ("Array size of PUT is too small.");
697 /* We copy the seed given by the user. */
698 for (i
= 0; i
< kiss_size
; i
++)
699 memcpy (seed
+ i
* sizeof(GFC_UINTEGER_4
),
700 &(put
->data
[(kiss_size
- 1 - i
) * GFC_DESCRIPTOR_STRIDE(put
,0)]),
701 sizeof(GFC_UINTEGER_4
));
703 /* We put it after scrambling the bytes, to paper around users who
704 provide seeds with quality only in the lower or upper part. */
705 scramble_seed ((unsigned char *) kiss_seed
, seed
, 4*kiss_size
);
708 /* Return the seed to GET data. */
711 /* If the rank of the array is not 1, abort. */
712 if (GFC_DESCRIPTOR_RANK (get
) != 1)
713 runtime_error ("Array rank of GET is not 1.");
715 /* If the array is too small, abort. */
716 if (GFC_DESCRIPTOR_EXTENT(get
,0) < kiss_size
)
717 runtime_error ("Array size of GET is too small.");
719 /* Unscramble the seed. */
720 unscramble_seed (seed
, (unsigned char *) kiss_seed
, 4*kiss_size
);
722 /* Then copy it back to the user variable. */
723 for (i
= 0; i
< kiss_size
; i
++)
724 memcpy (&(get
->data
[(kiss_size
- 1 - i
) * GFC_DESCRIPTOR_STRIDE(get
,0)]),
725 seed
+ i
* sizeof(GFC_UINTEGER_4
),
726 sizeof(GFC_UINTEGER_4
));
729 __gthread_mutex_unlock (&random_lock
);
731 iexport(random_seed_i4
);
735 random_seed_i8 (GFC_INTEGER_8
*size
, gfc_array_i8
*put
, gfc_array_i8
*get
)
739 __gthread_mutex_lock (&random_lock
);
741 /* Check that we only have one argument present. */
742 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
743 runtime_error ("RANDOM_SEED should have at most one argument present.");
745 /* From the standard: "If no argument is present, the processor assigns
746 a processor-dependent value to the seed." */
747 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
748 for (i
= 0; i
< kiss_size
; i
++)
749 kiss_seed
[i
] = kiss_default_seed
[i
];
752 *size
= kiss_size
/ 2;
756 /* If the rank of the array is not 1, abort. */
757 if (GFC_DESCRIPTOR_RANK (put
) != 1)
758 runtime_error ("Array rank of PUT is not 1.");
760 /* If the array is too small, abort. */
761 if (GFC_DESCRIPTOR_EXTENT(put
,0) < kiss_size
/ 2)
762 runtime_error ("Array size of PUT is too small.");
764 /* This code now should do correct strides. */
765 for (i
= 0; i
< kiss_size
/ 2; i
++)
766 memcpy (&kiss_seed
[2*i
], &(put
->data
[i
* GFC_DESCRIPTOR_STRIDE(put
,0)]),
767 sizeof (GFC_UINTEGER_8
));
770 /* Return the seed to GET data. */
773 /* If the rank of the array is not 1, abort. */
774 if (GFC_DESCRIPTOR_RANK (get
) != 1)
775 runtime_error ("Array rank of GET is not 1.");
777 /* If the array is too small, abort. */
778 if (GFC_DESCRIPTOR_EXTENT(get
,0) < kiss_size
/ 2)
779 runtime_error ("Array size of GET is too small.");
781 /* This code now should do correct strides. */
782 for (i
= 0; i
< kiss_size
/ 2; i
++)
783 memcpy (&(get
->data
[i
* GFC_DESCRIPTOR_STRIDE(get
,0)]), &kiss_seed
[2*i
],
784 sizeof (GFC_UINTEGER_8
));
787 __gthread_mutex_unlock (&random_lock
);
789 iexport(random_seed_i8
);
792 #ifndef __GTHREAD_MUTEX_INIT
793 static void __attribute__((constructor
))
796 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock
);