1 /* Implementation of the RANDOM intrinsics
2 Copyright 2002, 2004, 2005, 2006, 2007 Free Software Foundation, Inc.
3 Contributed by Lars Segerlund <seger@linuxmail.org>
6 This file is part of the GNU Fortran 95 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 2 of the License, or (at your option) any later version.
13 In addition to the permissions in the GNU General Public License, the
14 Free Software Foundation gives you unlimited permission to link the
15 compiled version of this file into combinations with other programs,
16 and to distribute those combinations without any restriction coming
17 from the use of this file. (The General Public License restrictions
18 do apply in other respects; for example, they cover modification of
19 the file, and distribution when not linked into a combine
22 Ligbfortran is distributed in the hope that it will be useful,
23 but WITHOUT ANY WARRANTY; without even the implied warranty of
24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 GNU General Public License for more details.
27 You should have received a copy of the GNU General Public
28 License along with libgfortran; see the file COPYING. If not,
29 write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
30 Boston, MA 02110-1301, USA. */
32 #include "libgfortran.h"
36 extern void random_r4 (GFC_REAL_4
*);
37 iexport_proto(random_r4
);
39 extern void random_r8 (GFC_REAL_8
*);
40 iexport_proto(random_r8
);
42 extern void arandom_r4 (gfc_array_r4
*);
43 export_proto(arandom_r4
);
45 extern void arandom_r8 (gfc_array_r8
*);
46 export_proto(arandom_r8
);
48 #ifdef HAVE_GFC_REAL_10
50 extern void random_r10 (GFC_REAL_10
*);
51 iexport_proto(random_r10
);
53 extern void arandom_r10 (gfc_array_r10
*);
54 export_proto(arandom_r10
);
58 #ifdef HAVE_GFC_REAL_16
60 extern void random_r16 (GFC_REAL_16
*);
61 iexport_proto(random_r16
);
63 extern void arandom_r16 (gfc_array_r16
*);
64 export_proto(arandom_r16
);
68 #ifdef __GTHREAD_MUTEX_INIT
69 static __gthread_mutex_t random_lock
= __GTHREAD_MUTEX_INIT
;
71 static __gthread_mutex_t random_lock
;
74 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
75 GFC_REAL_* types in the range of [0,1). If GFC_REAL_*_RADIX are 2
76 or 16, respectively, we mask off the bits that don't fit into the
77 correct GFC_REAL_*, convert to the real type, then multiply by the
83 rnumber_4 (GFC_REAL_4
*f
, GFC_UINTEGER_4 v
)
86 #if GFC_REAL_4_RADIX == 2
87 mask
= ~ (GFC_UINTEGER_4
) 0u << (32 - GFC_REAL_4_DIGITS
);
88 #elif GFC_REAL_4_RADIX == 16
89 mask
= ~ (GFC_UINTEGER_4
) 0u << ((8 - GFC_REAL_4_DIGITS
) * 4);
91 #error "GFC_REAL_4_RADIX has unknown value"
94 *f
= (GFC_REAL_4
) v
* (GFC_REAL_4
) 0x1.p
-32f
;
98 rnumber_8 (GFC_REAL_8
*f
, GFC_UINTEGER_8 v
)
101 #if GFC_REAL_8_RADIX == 2
102 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_8_DIGITS
);
103 #elif GFC_REAL_8_RADIX == 16
104 mask
= ~ (GFC_UINTEGER_8
) 0u << (16 - GFC_REAL_8_DIGITS
) * 4);
106 #error "GFC_REAL_8_RADIX has unknown value"
109 *f
= (GFC_REAL_8
) v
* (GFC_REAL_8
) 0x1.p
-64;
112 #ifdef HAVE_GFC_REAL_10
115 rnumber_10 (GFC_REAL_10
*f
, GFC_UINTEGER_8 v
)
118 #if GFC_REAL_10_RADIX == 2
119 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_10_DIGITS
);
120 #elif GFC_REAL_10_RADIX == 16
121 mask
= ~ (GFC_UINTEGER_10
) 0u << ((16 - GFC_REAL_10_DIGITS
) * 4);
123 #error "GFC_REAL_10_RADIX has unknown value"
126 *f
= (GFC_REAL_10
) v
* (GFC_REAL_10
) 0x1.p
-64;
130 #ifdef HAVE_GFC_REAL_16
132 /* For REAL(KIND=16), we only need to mask off the lower bits. */
135 rnumber_16 (GFC_REAL_16
*f
, GFC_UINTEGER_8 v1
, GFC_UINTEGER_8 v2
)
138 #if GFC_REAL_16_RADIX == 2
139 mask
= ~ (GFC_UINTEGER_8
) 0u << (128 - GFC_REAL_16_DIGITS
);
140 #elif GFC_REAL_16_RADIX == 16
141 mask
= ~ (GFC_UINTEGER_8
) 0u << ((32 - GFC_REAL_16_DIGITS
) * 4);
143 #error "GFC_REAL_16_RADIX has unknown value"
146 *f
= (GFC_REAL_16
) v1
* (GFC_REAL_16
) 0x1.p
-64
147 + (GFC_REAL_16
) v2
* (GFC_REAL_16
) 0x1.p
-128;
150 /* libgfortran previously had a Mersenne Twister, taken from the paper:
152 Mersenne Twister: 623-dimensionally equidistributed
153 uniform pseudorandom generator.
155 by Makoto Matsumoto & Takuji Nishimura
156 which appeared in the: ACM Transactions on Modelling and Computer
157 Simulations: Special Issue on Uniform Random Number
158 Generation. ( Early in 1998 ).
160 The Mersenne Twister code was replaced due to
162 (1) Simple user specified seeds lead to really bad sequences for
163 nearly 100000 random numbers.
164 (2) open(), read(), and close() were not properly declared via header
166 (3) The global index i was abused and caused unexpected behavior with
171 libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
172 random number generator. This PRNG combines:
174 (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
176 (2) A 3-shift shift-register generator with a period of 2^32-1,
177 (3) Two 16-bit multiply-with-carry generators with a period of
178 597273182964842497 > 2^59.
180 The overall period exceeds 2^123.
182 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
184 The above web site has an archive of a newsgroup posting from George
185 Marsaglia with the statement:
187 Subject: Random numbers for C: Improvements.
188 Date: Fri, 15 Jan 1999 11:41:47 -0500
189 From: George Marsaglia <geo@stat.fsu.edu>
190 Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
191 References: <369B5E30.65A55FD1@stat.fsu.edu>
192 Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
195 As I hoped, several suggestions have led to
196 improvements in the code for RNG's I proposed for
197 use in C. (See the thread "Random numbers for C: Some
198 suggestions" in previous postings.) The improved code
201 A question of copyright has also been raised. Unlike
202 DIEHARD, there is no copyright on the code below. You
203 are free to use it in any way you want, but you may
204 wish to acknowledge the source, as a courtesy.
206 "There is no copyright on the code below." included the original
209 /* We use three KISS random number generators, with different
211 As a matter of Quality of Implementation, the random numbers
212 we generate for different REAL kinds, starting from the same
213 seed, are always the same up to the precision of these types.
214 We do this by using three generators with different seeds, the
215 first one always for the most significant bits, the second one
216 for bits 33..64 (if present in the REAL kind), and the third one
217 (called twice) for REAL(16).
220 #define GFC_SL(k, n) ((k)^((k)<<(n)))
221 #define GFC_SR(k, n) ((k)^((k)>>(n)))
223 /* Reference for the seed:
224 From: "George Marsaglia" <g...@stat.fsu.edu>
226 Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
228 The KISS RNG uses four seeds, x, y, z, c,
229 with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
230 except that the two pairs
231 z=0,c=0 and z=2^32-1,c=698769068
235 #define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
236 #define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
237 #ifdef HAVE_GFC_REAL_16
238 #define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
241 static GFC_UINTEGER_4 kiss_seed
[] = {
244 #ifdef HAVE_GFC_REAL_16
249 static GFC_UINTEGER_4 kiss_default_seed
[] = {
252 #ifdef HAVE_GFC_REAL_16
257 static const GFC_INTEGER_4 kiss_size
= sizeof(kiss_seed
)/sizeof(kiss_seed
[0]);
259 static GFC_UINTEGER_4
* const kiss_seed_1
= kiss_seed
;
260 static GFC_UINTEGER_4
* const kiss_seed_2
= kiss_seed
+ 4;
262 #ifdef HAVE_GFC_REAL_16
263 static GFC_UINTEGER_4
* const kiss_seed_3
= kiss_seed
+ 8;
266 /* kiss_random_kernel() returns an integer value in the range of
267 (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
268 should be uniform. */
270 static GFC_UINTEGER_4
271 kiss_random_kernel(GFC_UINTEGER_4
* seed
)
275 seed
[0] = 69069 * seed
[0] + 1327217885;
276 seed
[1] = GFC_SL(GFC_SR(GFC_SL(seed
[1],13),17),5);
277 seed
[2] = 18000 * (seed
[2] & 65535) + (seed
[2] >> 16);
278 seed
[3] = 30903 * (seed
[3] & 65535) + (seed
[3] >> 16);
279 kiss
= seed
[0] + seed
[1] + (seed
[2] << 16) + seed
[3];
284 /* This function produces a REAL(4) value from the uniform distribution
288 random_r4 (GFC_REAL_4
*x
)
292 __gthread_mutex_lock (&random_lock
);
293 kiss
= kiss_random_kernel (kiss_seed_1
);
295 __gthread_mutex_unlock (&random_lock
);
299 /* This function produces a REAL(8) value from the uniform distribution
303 random_r8 (GFC_REAL_8
*x
)
307 __gthread_mutex_lock (&random_lock
);
308 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
309 kiss
+= kiss_random_kernel (kiss_seed_2
);
311 __gthread_mutex_unlock (&random_lock
);
315 #ifdef HAVE_GFC_REAL_10
317 /* This function produces a REAL(10) value from the uniform distribution
321 random_r10 (GFC_REAL_10
*x
)
325 __gthread_mutex_lock (&random_lock
);
326 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
327 kiss
+= kiss_random_kernel (kiss_seed_2
);
328 rnumber_10 (x
, kiss
);
329 __gthread_mutex_unlock (&random_lock
);
335 /* This function produces a REAL(16) value from the uniform distribution
338 #ifdef HAVE_GFC_REAL_16
341 random_r16 (GFC_REAL_16
*x
)
343 GFC_UINTEGER_8 kiss1
, kiss2
;
345 __gthread_mutex_lock (&random_lock
);
346 kiss1
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
347 kiss1
+= kiss_random_kernel (kiss_seed_2
);
349 kiss2
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_3
)) << 32;
350 kiss2
+= kiss_random_kernel (kiss_seed_3
);
352 rnumber_16 (x
, kiss1
, kiss2
);
353 __gthread_mutex_unlock (&random_lock
);
359 /* This function fills a REAL(4) array with values from the uniform
360 distribution with range [0,1). */
363 arandom_r4 (gfc_array_r4
*x
)
365 index_type count
[GFC_MAX_DIMENSIONS
];
366 index_type extent
[GFC_MAX_DIMENSIONS
];
367 index_type stride
[GFC_MAX_DIMENSIONS
];
376 dim
= GFC_DESCRIPTOR_RANK (x
);
378 for (n
= 0; n
< dim
; n
++)
381 stride
[n
] = x
->dim
[n
].stride
;
382 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
389 __gthread_mutex_lock (&random_lock
);
393 /* random_r4 (dest); */
394 kiss
= kiss_random_kernel (kiss_seed_1
);
395 rnumber_4 (dest
, kiss
);
397 /* Advance to the next element. */
400 /* Advance to the next source element. */
402 while (count
[n
] == extent
[n
])
404 /* When we get to the end of a dimension, reset it and increment
405 the next dimension. */
407 /* We could precalculate these products, but this is a less
408 frequently used path so probably not worth it. */
409 dest
-= stride
[n
] * extent
[n
];
423 __gthread_mutex_unlock (&random_lock
);
426 /* This function fills a REAL(8) array with values from the uniform
427 distribution with range [0,1). */
430 arandom_r8 (gfc_array_r8
*x
)
432 index_type count
[GFC_MAX_DIMENSIONS
];
433 index_type extent
[GFC_MAX_DIMENSIONS
];
434 index_type stride
[GFC_MAX_DIMENSIONS
];
443 dim
= GFC_DESCRIPTOR_RANK (x
);
445 for (n
= 0; n
< dim
; n
++)
448 stride
[n
] = x
->dim
[n
].stride
;
449 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
456 __gthread_mutex_lock (&random_lock
);
460 /* random_r8 (dest); */
461 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
462 kiss
+= kiss_random_kernel (kiss_seed_2
);
463 rnumber_8 (dest
, kiss
);
465 /* Advance to the next element. */
468 /* Advance to the next source element. */
470 while (count
[n
] == extent
[n
])
472 /* When we get to the end of a dimension, reset it and increment
473 the next dimension. */
475 /* We could precalculate these products, but this is a less
476 frequently used path so probably not worth it. */
477 dest
-= stride
[n
] * extent
[n
];
491 __gthread_mutex_unlock (&random_lock
);
494 #ifdef HAVE_GFC_REAL_10
496 /* This function fills a REAL(10) array with values from the uniform
497 distribution with range [0,1). */
500 arandom_r10 (gfc_array_r10
*x
)
502 index_type count
[GFC_MAX_DIMENSIONS
];
503 index_type extent
[GFC_MAX_DIMENSIONS
];
504 index_type stride
[GFC_MAX_DIMENSIONS
];
513 dim
= GFC_DESCRIPTOR_RANK (x
);
515 for (n
= 0; n
< dim
; n
++)
518 stride
[n
] = x
->dim
[n
].stride
;
519 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
526 __gthread_mutex_lock (&random_lock
);
530 /* random_r10 (dest); */
531 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
532 kiss
+= kiss_random_kernel (kiss_seed_2
);
533 rnumber_10 (dest
, kiss
);
535 /* Advance to the next element. */
538 /* Advance to the next source element. */
540 while (count
[n
] == extent
[n
])
542 /* When we get to the end of a dimension, reset it and increment
543 the next dimension. */
545 /* We could precalculate these products, but this is a less
546 frequently used path so probably not worth it. */
547 dest
-= stride
[n
] * extent
[n
];
561 __gthread_mutex_unlock (&random_lock
);
566 #ifdef HAVE_GFC_REAL_16
568 /* This function fills a REAL(16) array with values from the uniform
569 distribution with range [0,1). */
572 arandom_r16 (gfc_array_r16
*x
)
574 index_type count
[GFC_MAX_DIMENSIONS
];
575 index_type extent
[GFC_MAX_DIMENSIONS
];
576 index_type stride
[GFC_MAX_DIMENSIONS
];
580 GFC_UINTEGER_8 kiss1
, kiss2
;
585 dim
= GFC_DESCRIPTOR_RANK (x
);
587 for (n
= 0; n
< dim
; n
++)
590 stride
[n
] = x
->dim
[n
].stride
;
591 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
598 __gthread_mutex_lock (&random_lock
);
602 /* random_r16 (dest); */
603 kiss1
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
604 kiss1
+= kiss_random_kernel (kiss_seed_2
);
606 kiss2
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_3
)) << 32;
607 kiss2
+= kiss_random_kernel (kiss_seed_3
);
609 rnumber_16 (dest
, kiss1
, kiss2
);
611 /* Advance to the next element. */
614 /* Advance to the next source element. */
616 while (count
[n
] == extent
[n
])
618 /* When we get to the end of a dimension, reset it and increment
619 the next dimension. */
621 /* We could precalculate these products, but this is a less
622 frequently used path so probably not worth it. */
623 dest
-= stride
[n
] * extent
[n
];
637 __gthread_mutex_unlock (&random_lock
);
645 scramble_seed (unsigned char *dest
, unsigned char *src
, int size
)
649 for (i
= 0; i
< size
; i
++)
650 dest
[(i
% 2) * (size
/ 2) + i
/ 2] = src
[i
];
655 unscramble_seed (unsigned char *dest
, unsigned char *src
, int size
)
659 for (i
= 0; i
< size
; i
++)
660 dest
[i
] = src
[(i
% 2) * (size
/ 2) + i
/ 2];
665 /* random_seed is used to seed the PRNG with either a default
666 set of seeds or user specified set of seeds. random_seed
667 must be called with no argument or exactly one argument. */
670 random_seed_i4 (GFC_INTEGER_4
*size
, gfc_array_i4
*put
, gfc_array_i4
*get
)
673 unsigned char seed
[4*kiss_size
];
675 __gthread_mutex_lock (&random_lock
);
677 /* Check that we only have one argument present. */
678 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
679 runtime_error ("RANDOM_SEED should have at most one argument present.");
681 /* From the standard: "If no argument is present, the processor assigns
682 a processor-dependent value to the seed." */
683 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
684 for (i
= 0; i
< kiss_size
; i
++)
685 kiss_seed
[i
] = kiss_default_seed
[i
];
692 /* If the rank of the array is not 1, abort. */
693 if (GFC_DESCRIPTOR_RANK (put
) != 1)
694 runtime_error ("Array rank of PUT is not 1.");
696 /* If the array is too small, abort. */
697 if (((put
->dim
[0].ubound
+ 1 - put
->dim
[0].lbound
)) < kiss_size
)
698 runtime_error ("Array size of PUT is too small.");
700 /* We copy the seed given by the user. */
701 for (i
= 0; i
< kiss_size
; i
++)
702 memcpy (seed
+ i
* sizeof(GFC_UINTEGER_4
),
703 &(put
->data
[(kiss_size
- 1 - i
) * put
->dim
[0].stride
]),
704 sizeof(GFC_UINTEGER_4
));
706 /* We put it after scrambling the bytes, to paper around users who
707 provide seeds with quality only in the lower or upper part. */
708 scramble_seed ((unsigned char *) kiss_seed
, seed
, 4*kiss_size
);
711 /* Return the seed to GET data. */
714 /* If the rank of the array is not 1, abort. */
715 if (GFC_DESCRIPTOR_RANK (get
) != 1)
716 runtime_error ("Array rank of GET is not 1.");
718 /* If the array is too small, abort. */
719 if (((get
->dim
[0].ubound
+ 1 - get
->dim
[0].lbound
)) < kiss_size
)
720 runtime_error ("Array size of GET is too small.");
722 /* Unscramble the seed. */
723 unscramble_seed (seed
, (unsigned char *) kiss_seed
, 4*kiss_size
);
725 /* Then copy it back to the user variable. */
726 for (i
= 0; i
< kiss_size
; i
++)
727 memcpy (&(get
->data
[(kiss_size
- 1 - i
) * get
->dim
[0].stride
]),
728 seed
+ i
* sizeof(GFC_UINTEGER_4
),
729 sizeof(GFC_UINTEGER_4
));
732 __gthread_mutex_unlock (&random_lock
);
734 iexport(random_seed_i4
);
738 random_seed_i8 (GFC_INTEGER_8
*size
, gfc_array_i8
*put
, gfc_array_i8
*get
)
742 __gthread_mutex_lock (&random_lock
);
744 /* Check that we only have one argument present. */
745 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
746 runtime_error ("RANDOM_SEED should have at most one argument present.");
748 /* From the standard: "If no argument is present, the processor assigns
749 a processor-dependent value to the seed." */
750 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
751 for (i
= 0; i
< kiss_size
; i
++)
752 kiss_seed
[i
] = kiss_default_seed
[i
];
755 *size
= kiss_size
/ 2;
759 /* If the rank of the array is not 1, abort. */
760 if (GFC_DESCRIPTOR_RANK (put
) != 1)
761 runtime_error ("Array rank of PUT is not 1.");
763 /* If the array is too small, abort. */
764 if (((put
->dim
[0].ubound
+ 1 - put
->dim
[0].lbound
)) < kiss_size
/ 2)
765 runtime_error ("Array size of PUT is too small.");
767 /* This code now should do correct strides. */
768 for (i
= 0; i
< kiss_size
/ 2; i
++)
769 memcpy (&kiss_seed
[2*i
], &(put
->data
[i
* put
->dim
[0].stride
]),
770 sizeof (GFC_UINTEGER_8
));
773 /* Return the seed to GET data. */
776 /* If the rank of the array is not 1, abort. */
777 if (GFC_DESCRIPTOR_RANK (get
) != 1)
778 runtime_error ("Array rank of GET is not 1.");
780 /* If the array is too small, abort. */
781 if (((get
->dim
[0].ubound
+ 1 - get
->dim
[0].lbound
)) < kiss_size
/ 2)
782 runtime_error ("Array size of GET is too small.");
784 /* This code now should do correct strides. */
785 for (i
= 0; i
< kiss_size
/ 2; i
++)
786 memcpy (&(get
->data
[i
* get
->dim
[0].stride
]), &kiss_seed
[2*i
],
787 sizeof (GFC_UINTEGER_8
));
790 __gthread_mutex_unlock (&random_lock
);
792 iexport(random_seed_i8
);
795 #ifndef __GTHREAD_MUTEX_INIT
796 static void __attribute__((constructor
))
799 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock
);