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. */
33 #include "libgfortran.h"
37 extern void random_r4 (GFC_REAL_4
*);
38 iexport_proto(random_r4
);
40 extern void random_r8 (GFC_REAL_8
*);
41 iexport_proto(random_r8
);
43 extern void arandom_r4 (gfc_array_r4
*);
44 export_proto(arandom_r4
);
46 extern void arandom_r8 (gfc_array_r8
*);
47 export_proto(arandom_r8
);
49 #ifdef HAVE_GFC_REAL_10
51 extern void random_r10 (GFC_REAL_10
*);
52 iexport_proto(random_r10
);
54 extern void arandom_r10 (gfc_array_r10
*);
55 export_proto(arandom_r10
);
59 #ifdef HAVE_GFC_REAL_16
61 extern void random_r16 (GFC_REAL_16
*);
62 iexport_proto(random_r16
);
64 extern void arandom_r16 (gfc_array_r16
*);
65 export_proto(arandom_r16
);
69 #ifdef __GTHREAD_MUTEX_INIT
70 static __gthread_mutex_t random_lock
= __GTHREAD_MUTEX_INIT
;
72 static __gthread_mutex_t random_lock
;
75 /* Helper routines to map a GFC_UINTEGER_* to the corresponding
76 GFC_REAL_* types in the range of [0,1). If GFC_REAL_*_RADIX are 2
77 or 16, respectively, we mask off the bits that don't fit into the
78 correct GFC_REAL_*, convert to the real type, then multiply by the
84 rnumber_4 (GFC_REAL_4
*f
, GFC_UINTEGER_4 v
)
87 #if GFC_REAL_4_RADIX == 2
88 mask
= ~ (GFC_UINTEGER_4
) 0u << (32 - GFC_REAL_4_DIGITS
);
89 #elif GFC_REAL_4_RADIX == 16
90 mask
= ~ (GFC_UINTEGER_4
) 0u << ((8 - GFC_REAL_4_DIGITS
) * 4);
92 #error "GFC_REAL_4_RADIX has unknown value"
95 *f
= (GFC_REAL_4
) v
* (GFC_REAL_4
) 0x1.p
-32f
;
99 rnumber_8 (GFC_REAL_8
*f
, GFC_UINTEGER_8 v
)
102 #if GFC_REAL_8_RADIX == 2
103 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_8_DIGITS
);
104 #elif GFC_REAL_8_RADIX == 16
105 mask
= ~ (GFC_UINTEGER_8
) 0u << (16 - GFC_REAL_8_DIGITS
) * 4);
107 #error "GFC_REAL_8_RADIX has unknown value"
110 *f
= (GFC_REAL_8
) v
* (GFC_REAL_8
) 0x1.p
-64;
113 #ifdef HAVE_GFC_REAL_10
116 rnumber_10 (GFC_REAL_10
*f
, GFC_UINTEGER_8 v
)
119 #if GFC_REAL_10_RADIX == 2
120 mask
= ~ (GFC_UINTEGER_8
) 0u << (64 - GFC_REAL_10_DIGITS
);
121 #elif GFC_REAL_10_RADIX == 16
122 mask
= ~ (GFC_UINTEGER_10
) 0u << ((16 - GFC_REAL_10_DIGITS
) * 4);
124 #error "GFC_REAL_10_RADIX has unknown value"
127 *f
= (GFC_REAL_10
) v
* (GFC_REAL_10
) 0x1.p
-64;
131 #ifdef HAVE_GFC_REAL_16
133 /* For REAL(KIND=16), we only need to mask off the lower bits. */
136 rnumber_16 (GFC_REAL_16
*f
, GFC_UINTEGER_8 v1
, GFC_UINTEGER_8 v2
)
139 #if GFC_REAL_16_RADIX == 2
140 mask
= ~ (GFC_UINTEGER_8
) 0u << (128 - GFC_REAL_16_DIGITS
);
141 #elif GFC_REAL_16_RADIX == 16
142 mask
= ~ (GFC_UINTEGER_8
) 0u << ((32 - GFC_REAL_16_DIGITS
) * 4);
144 #error "GFC_REAL_16_RADIX has unknown value"
147 *f
= (GFC_REAL_16
) v1
* (GFC_REAL_16
) 0x1.p
-64
148 + (GFC_REAL_16
) v2
* (GFC_REAL_16
) 0x1.p
-128;
151 /* libgfortran previously had a Mersenne Twister, taken from the paper:
153 Mersenne Twister: 623-dimensionally equidistributed
154 uniform pseudorandom generator.
156 by Makoto Matsumoto & Takuji Nishimura
157 which appeared in the: ACM Transactions on Modelling and Computer
158 Simulations: Special Issue on Uniform Random Number
159 Generation. ( Early in 1998 ).
161 The Mersenne Twister code was replaced due to
163 (1) Simple user specified seeds lead to really bad sequences for
164 nearly 100000 random numbers.
165 (2) open(), read(), and close() were not properly declared via header
167 (3) The global index i was abused and caused unexpected behavior with
172 libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
173 random number generator. This PRNG combines:
175 (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
177 (2) A 3-shift shift-register generator with a period of 2^32-1,
178 (3) Two 16-bit multiply-with-carry generators with a period of
179 597273182964842497 > 2^59.
181 The overall period exceeds 2^123.
183 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
185 The above web site has an archive of a newsgroup posting from George
186 Marsaglia with the statement:
188 Subject: Random numbers for C: Improvements.
189 Date: Fri, 15 Jan 1999 11:41:47 -0500
190 From: George Marsaglia <geo@stat.fsu.edu>
191 Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
192 References: <369B5E30.65A55FD1@stat.fsu.edu>
193 Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
196 As I hoped, several suggestions have led to
197 improvements in the code for RNG's I proposed for
198 use in C. (See the thread "Random numbers for C: Some
199 suggestions" in previous postings.) The improved code
202 A question of copyright has also been raised. Unlike
203 DIEHARD, there is no copyright on the code below. You
204 are free to use it in any way you want, but you may
205 wish to acknowledge the source, as a courtesy.
207 "There is no copyright on the code below." included the original
210 /* We use three KISS random number generators, with different
212 As a matter of Quality of Implementation, the random numbers
213 we generate for different REAL kinds, starting from the same
214 seed, are always the same up to the precision of these types.
215 We do this by using three generators with different seeds, the
216 first one always for the most significant bits, the second one
217 for bits 33..64 (if present in the REAL kind), and the third one
218 (called twice) for REAL(16).
221 #define GFC_SL(k, n) ((k)^((k)<<(n)))
222 #define GFC_SR(k, n) ((k)^((k)>>(n)))
224 /* Reference for the seed:
225 From: "George Marsaglia" <g...@stat.fsu.edu>
227 Message-ID: <e7CcnWxczriWssCjXTWc3A@comcast.com>
229 The KISS RNG uses four seeds, x, y, z, c,
230 with 0<=x<2^32, 0<y<2^32, 0<=z<2^32, 0<=c<698769069
231 except that the two pairs
232 z=0,c=0 and z=2^32-1,c=698769068
236 #define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069
237 #define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021
238 #ifdef HAVE_GFC_REAL_16
239 #define KISS_DEFAULT_SEED_3 573658661, 185639104, 582619469, 296736107
242 static GFC_UINTEGER_4 kiss_seed
[] = {
245 #ifdef HAVE_GFC_REAL_16
250 static GFC_UINTEGER_4 kiss_default_seed
[] = {
253 #ifdef HAVE_GFC_REAL_16
258 static const GFC_INTEGER_4 kiss_size
= sizeof(kiss_seed
)/sizeof(kiss_seed
[0]);
260 static GFC_UINTEGER_4
* const kiss_seed_1
= kiss_seed
;
261 static GFC_UINTEGER_4
* const kiss_seed_2
= kiss_seed
+ 4;
263 #ifdef HAVE_GFC_REAL_16
264 static GFC_UINTEGER_4
* const kiss_seed_3
= kiss_seed
+ 8;
267 /* kiss_random_kernel() returns an integer value in the range of
268 (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
269 should be uniform. */
271 static GFC_UINTEGER_4
272 kiss_random_kernel(GFC_UINTEGER_4
* seed
)
276 seed
[0] = 69069 * seed
[0] + 1327217885;
277 seed
[1] = GFC_SL(GFC_SR(GFC_SL(seed
[1],13),17),5);
278 seed
[2] = 18000 * (seed
[2] & 65535) + (seed
[2] >> 16);
279 seed
[3] = 30903 * (seed
[3] & 65535) + (seed
[3] >> 16);
280 kiss
= seed
[0] + seed
[1] + (seed
[2] << 16) + seed
[3];
285 /* This function produces a REAL(4) value from the uniform distribution
289 random_r4 (GFC_REAL_4
*x
)
293 __gthread_mutex_lock (&random_lock
);
294 kiss
= kiss_random_kernel (kiss_seed_1
);
296 __gthread_mutex_unlock (&random_lock
);
300 /* This function produces a REAL(8) value from the uniform distribution
304 random_r8 (GFC_REAL_8
*x
)
308 __gthread_mutex_lock (&random_lock
);
309 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
310 kiss
+= kiss_random_kernel (kiss_seed_2
);
312 __gthread_mutex_unlock (&random_lock
);
316 #ifdef HAVE_GFC_REAL_10
318 /* This function produces a REAL(10) value from the uniform distribution
322 random_r10 (GFC_REAL_10
*x
)
326 __gthread_mutex_lock (&random_lock
);
327 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
328 kiss
+= kiss_random_kernel (kiss_seed_2
);
329 rnumber_10 (x
, kiss
);
330 __gthread_mutex_unlock (&random_lock
);
336 /* This function produces a REAL(16) value from the uniform distribution
339 #ifdef HAVE_GFC_REAL_16
342 random_r16 (GFC_REAL_16
*x
)
344 GFC_UINTEGER_8 kiss1
, kiss2
;
346 __gthread_mutex_lock (&random_lock
);
347 kiss1
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
348 kiss1
+= kiss_random_kernel (kiss_seed_2
);
350 kiss2
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_3
)) << 32;
351 kiss2
+= kiss_random_kernel (kiss_seed_3
);
353 rnumber_16 (x
, kiss1
, kiss2
);
354 __gthread_mutex_unlock (&random_lock
);
360 /* This function fills a REAL(4) array with values from the uniform
361 distribution with range [0,1). */
364 arandom_r4 (gfc_array_r4
*x
)
366 index_type count
[GFC_MAX_DIMENSIONS
];
367 index_type extent
[GFC_MAX_DIMENSIONS
];
368 index_type stride
[GFC_MAX_DIMENSIONS
];
377 dim
= GFC_DESCRIPTOR_RANK (x
);
379 for (n
= 0; n
< dim
; n
++)
382 stride
[n
] = x
->dim
[n
].stride
;
383 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
390 __gthread_mutex_lock (&random_lock
);
394 /* random_r4 (dest); */
395 kiss
= kiss_random_kernel (kiss_seed_1
);
396 rnumber_4 (dest
, kiss
);
398 /* Advance to the next element. */
401 /* Advance to the next source element. */
403 while (count
[n
] == extent
[n
])
405 /* When we get to the end of a dimension, reset it and increment
406 the next dimension. */
408 /* We could precalculate these products, but this is a less
409 frequently used path so probably not worth it. */
410 dest
-= stride
[n
] * extent
[n
];
424 __gthread_mutex_unlock (&random_lock
);
427 /* This function fills a REAL(8) array with values from the uniform
428 distribution with range [0,1). */
431 arandom_r8 (gfc_array_r8
*x
)
433 index_type count
[GFC_MAX_DIMENSIONS
];
434 index_type extent
[GFC_MAX_DIMENSIONS
];
435 index_type stride
[GFC_MAX_DIMENSIONS
];
444 dim
= GFC_DESCRIPTOR_RANK (x
);
446 for (n
= 0; n
< dim
; n
++)
449 stride
[n
] = x
->dim
[n
].stride
;
450 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
457 __gthread_mutex_lock (&random_lock
);
461 /* random_r8 (dest); */
462 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
463 kiss
+= kiss_random_kernel (kiss_seed_2
);
464 rnumber_8 (dest
, kiss
);
466 /* Advance to the next element. */
469 /* Advance to the next source element. */
471 while (count
[n
] == extent
[n
])
473 /* When we get to the end of a dimension, reset it and increment
474 the next dimension. */
476 /* We could precalculate these products, but this is a less
477 frequently used path so probably not worth it. */
478 dest
-= stride
[n
] * extent
[n
];
492 __gthread_mutex_unlock (&random_lock
);
495 #ifdef HAVE_GFC_REAL_10
497 /* This function fills a REAL(10) array with values from the uniform
498 distribution with range [0,1). */
501 arandom_r10 (gfc_array_r10
*x
)
503 index_type count
[GFC_MAX_DIMENSIONS
];
504 index_type extent
[GFC_MAX_DIMENSIONS
];
505 index_type stride
[GFC_MAX_DIMENSIONS
];
514 dim
= GFC_DESCRIPTOR_RANK (x
);
516 for (n
= 0; n
< dim
; n
++)
519 stride
[n
] = x
->dim
[n
].stride
;
520 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
527 __gthread_mutex_lock (&random_lock
);
531 /* random_r10 (dest); */
532 kiss
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
533 kiss
+= kiss_random_kernel (kiss_seed_2
);
534 rnumber_10 (dest
, kiss
);
536 /* Advance to the next element. */
539 /* Advance to the next source element. */
541 while (count
[n
] == extent
[n
])
543 /* When we get to the end of a dimension, reset it and increment
544 the next dimension. */
546 /* We could precalculate these products, but this is a less
547 frequently used path so probably not worth it. */
548 dest
-= stride
[n
] * extent
[n
];
562 __gthread_mutex_unlock (&random_lock
);
567 #ifdef HAVE_GFC_REAL_16
569 /* This function fills a REAL(16) array with values from the uniform
570 distribution with range [0,1). */
573 arandom_r16 (gfc_array_r16
*x
)
575 index_type count
[GFC_MAX_DIMENSIONS
];
576 index_type extent
[GFC_MAX_DIMENSIONS
];
577 index_type stride
[GFC_MAX_DIMENSIONS
];
581 GFC_UINTEGER_8 kiss1
, kiss2
;
586 dim
= GFC_DESCRIPTOR_RANK (x
);
588 for (n
= 0; n
< dim
; n
++)
591 stride
[n
] = x
->dim
[n
].stride
;
592 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
599 __gthread_mutex_lock (&random_lock
);
603 /* random_r16 (dest); */
604 kiss1
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_1
)) << 32;
605 kiss1
+= kiss_random_kernel (kiss_seed_2
);
607 kiss2
= ((GFC_UINTEGER_8
) kiss_random_kernel (kiss_seed_3
)) << 32;
608 kiss2
+= kiss_random_kernel (kiss_seed_3
);
610 rnumber_16 (dest
, kiss1
, kiss2
);
612 /* Advance to the next element. */
615 /* Advance to the next source element. */
617 while (count
[n
] == extent
[n
])
619 /* When we get to the end of a dimension, reset it and increment
620 the next dimension. */
622 /* We could precalculate these products, but this is a less
623 frequently used path so probably not worth it. */
624 dest
-= stride
[n
] * extent
[n
];
638 __gthread_mutex_unlock (&random_lock
);
643 /* random_seed is used to seed the PRNG with either a default
644 set of seeds or user specified set of seeds. random_seed
645 must be called with no argument or exactly one argument. */
648 random_seed_i4 (GFC_INTEGER_4
*size
, gfc_array_i4
*put
, gfc_array_i4
*get
)
652 __gthread_mutex_lock (&random_lock
);
654 /* Check that we only have one argument present. */
655 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
656 runtime_error ("RANDOM_SEED should have at most one argument present.");
658 /* From the standard: "If no argument is present, the processor assigns
659 a processor-dependent value to the seed." */
660 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
661 for (i
= 0; i
< kiss_size
; i
++)
662 kiss_seed
[i
] = kiss_default_seed
[i
];
669 /* If the rank of the array is not 1, abort. */
670 if (GFC_DESCRIPTOR_RANK (put
) != 1)
671 runtime_error ("Array rank of PUT is not 1.");
673 /* If the array is too small, abort. */
674 if (((put
->dim
[0].ubound
+ 1 - put
->dim
[0].lbound
)) < kiss_size
)
675 runtime_error ("Array size of PUT is too small.");
677 /* This code now should do correct strides. */
678 for (i
= 0; i
< kiss_size
; i
++)
679 kiss_seed
[i
] = (GFC_UINTEGER_4
) put
->data
[i
* put
->dim
[0].stride
];
682 /* Return the seed to GET data. */
685 /* If the rank of the array is not 1, abort. */
686 if (GFC_DESCRIPTOR_RANK (get
) != 1)
687 runtime_error ("Array rank of GET is not 1.");
689 /* If the array is too small, abort. */
690 if (((get
->dim
[0].ubound
+ 1 - get
->dim
[0].lbound
)) < kiss_size
)
691 runtime_error ("Array size of GET is too small.");
693 /* This code now should do correct strides. */
694 for (i
= 0; i
< kiss_size
; i
++)
695 get
->data
[i
* get
->dim
[0].stride
] = (GFC_INTEGER_4
) kiss_seed
[i
];
698 __gthread_mutex_unlock (&random_lock
);
700 iexport(random_seed_i4
);
704 random_seed_i8 (GFC_INTEGER_8
*size
, gfc_array_i8
*put
, gfc_array_i8
*get
)
708 __gthread_mutex_lock (&random_lock
);
710 /* Check that we only have one argument present. */
711 if ((size
? 1 : 0) + (put
? 1 : 0) + (get
? 1 : 0) > 1)
712 runtime_error ("RANDOM_SEED should have at most one argument present.");
714 /* From the standard: "If no argument is present, the processor assigns
715 a processor-dependent value to the seed." */
716 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
717 for (i
= 0; i
< kiss_size
; i
++)
718 kiss_seed
[i
] = kiss_default_seed
[i
];
721 *size
= kiss_size
/ 2;
725 /* If the rank of the array is not 1, abort. */
726 if (GFC_DESCRIPTOR_RANK (put
) != 1)
727 runtime_error ("Array rank of PUT is not 1.");
729 /* If the array is too small, abort. */
730 if (((put
->dim
[0].ubound
+ 1 - put
->dim
[0].lbound
)) < kiss_size
/ 2)
731 runtime_error ("Array size of PUT is too small.");
733 /* This code now should do correct strides. */
734 for (i
= 0; i
< kiss_size
; i
+= 2)
735 memcpy (&kiss_seed
[i
], &(put
->data
[i
* put
->dim
[0].stride
]),
736 sizeof (GFC_UINTEGER_8
));
739 /* Return the seed to GET data. */
742 /* If the rank of the array is not 1, abort. */
743 if (GFC_DESCRIPTOR_RANK (get
) != 1)
744 runtime_error ("Array rank of GET is not 1.");
746 /* If the array is too small, abort. */
747 if (((get
->dim
[0].ubound
+ 1 - get
->dim
[0].lbound
)) < kiss_size
/ 2)
748 runtime_error ("Array size of GET is too small.");
750 /* This code now should do correct strides. */
751 for (i
= 0; i
< kiss_size
; i
+= 2)
752 memcpy (&(get
->data
[i
* get
->dim
[0].stride
]), &kiss_seed
[i
],
753 sizeof (GFC_UINTEGER_8
));
756 __gthread_mutex_unlock (&random_lock
);
758 iexport(random_seed_i8
);
761 #ifndef __GTHREAD_MUTEX_INIT
762 static void __attribute__((constructor
))
765 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock
);