1 /* Implementation of the RANDOM intrinsics
2 Copyright 2002, 2004, 2005 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"
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 __GTHREAD_MUTEX_INIT
49 static __gthread_mutex_t random_lock
= __GTHREAD_MUTEX_INIT
;
51 static __gthread_mutex_t random_lock
;
55 /* libgfortran previously had a Mersenne Twister, taken from the paper:
57 Mersenne Twister: 623-dimensionally equidistributed
58 uniform pseudorandom generator.
60 by Makoto Matsumoto & Takuji Nishimura
61 which appeared in the: ACM Transactions on Modelling and Computer
62 Simulations: Special Issue on Uniform Random Number
63 Generation. ( Early in 1998 ).
65 The Mersenne Twister code was replaced due to
67 (1) Simple user specified seeds lead to really bad sequences for
68 nearly 100000 random numbers.
69 (2) open(), read(), and close() were not properly declared via header
71 (3) The global index i was abused and caused unexpected behavior with
76 libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid)
77 random number generator. This PRNG combines:
79 (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
81 (2) A 3-shift shift-register generator with a period of 2^32-1,
82 (3) Two 16-bit multiply-with-carry generators with a period of
83 597273182964842497 > 2^59.
85 The overall period exceeds 2^123.
87 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
89 The above web site has an archive of a newsgroup posting from George
90 Marsaglia with the statement:
92 Subject: Random numbers for C: Improvements.
93 Date: Fri, 15 Jan 1999 11:41:47 -0500
94 From: George Marsaglia <geo@stat.fsu.edu>
95 Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
96 References: <369B5E30.65A55FD1@stat.fsu.edu>
97 Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
100 As I hoped, several suggestions have led to
101 improvements in the code for RNG's I proposed for
102 use in C. (See the thread "Random numbers for C: Some
103 suggestions" in previous postings.) The improved code
106 A question of copyright has also been raised. Unlike
107 DIEHARD, there is no copyright on the code below. You
108 are free to use it in any way you want, but you may
109 wish to acknowledge the source, as a courtesy.
111 "There is no copyright on the code below." included the original
114 #define GFC_SL(k, n) ((k)^((k)<<(n)))
115 #define GFC_SR(k, n) ((k)^((k)>>(n)))
117 static const GFC_INTEGER_4 kiss_size
= 4;
118 #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069}
119 static const GFC_UINTEGER_4 kiss_default_seed
[4] = KISS_DEFAULT_SEED
;
120 static GFC_UINTEGER_4 kiss_seed
[4] = KISS_DEFAULT_SEED
;
122 /* kiss_random_kernel() returns an integer value in the range of
123 (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
124 should be uniform. */
126 static GFC_UINTEGER_4
127 kiss_random_kernel(void)
131 kiss_seed
[0] = 69069 * kiss_seed
[0] + 1327217885;
132 kiss_seed
[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed
[1],13),17),5);
133 kiss_seed
[2] = 18000 * (kiss_seed
[2] & 65535) + (kiss_seed
[2] >> 16);
134 kiss_seed
[3] = 30903 * (kiss_seed
[3] & 65535) + (kiss_seed
[3] >> 16);
135 kiss
= kiss_seed
[0] + kiss_seed
[1] + (kiss_seed
[2] << 16) + kiss_seed
[3];
140 /* This function produces a REAL(4) value from the uniform distribution
144 random_r4 (GFC_REAL_4
*x
)
148 __gthread_mutex_lock (&random_lock
);
149 kiss
= kiss_random_kernel ();
150 /* Burn a random number, so the REAL*4 and REAL*8 functions
151 produce similar sequences of random numbers. */
152 kiss_random_kernel ();
153 *x
= normalize_r4_i4 (kiss
, ~(GFC_UINTEGER_4
) 0);
154 __gthread_mutex_unlock (&random_lock
);
158 /* This function produces a REAL(8) value from the uniform distribution
162 random_r8 (GFC_REAL_8
*x
)
166 __gthread_mutex_lock (&random_lock
);
167 kiss
= ((GFC_UINTEGER_8
)kiss_random_kernel ()) << 32;
168 kiss
+= kiss_random_kernel ();
169 *x
= normalize_r8_i8 (kiss
, ~(GFC_UINTEGER_8
) 0);
170 __gthread_mutex_unlock (&random_lock
);
174 /* This function fills a REAL(4) array with values from the uniform
175 distribution with range [0,1). */
178 arandom_r4 (gfc_array_r4
*x
)
180 index_type count
[GFC_MAX_DIMENSIONS
];
181 index_type extent
[GFC_MAX_DIMENSIONS
];
182 index_type stride
[GFC_MAX_DIMENSIONS
];
191 if (x
->dim
[0].stride
== 0)
192 x
->dim
[0].stride
= 1;
194 dim
= GFC_DESCRIPTOR_RANK (x
);
196 for (n
= 0; n
< dim
; n
++)
199 stride
[n
] = x
->dim
[n
].stride
;
200 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
207 __gthread_mutex_lock (&random_lock
);
211 /* random_r4 (dest); */
212 kiss
= kiss_random_kernel ();
213 /* Burn a random number, so the REAL*4 and REAL*8 functions
214 produce similar sequences of random numbers. */
215 kiss_random_kernel ();
216 *dest
= normalize_r4_i4 (kiss
, ~(GFC_UINTEGER_4
) 0);
218 /* Advance to the next element. */
221 /* Advance to the next source element. */
223 while (count
[n
] == extent
[n
])
225 /* When we get to the end of a dimension, reset it and increment
226 the next dimension. */
228 /* We could precalculate these products, but this is a less
229 frequently used path so probably not worth it. */
230 dest
-= stride
[n
] * extent
[n
];
244 __gthread_mutex_unlock (&random_lock
);
247 /* This function fills a REAL(8) array with values from the uniform
248 distribution with range [0,1). */
251 arandom_r8 (gfc_array_r8
*x
)
253 index_type count
[GFC_MAX_DIMENSIONS
];
254 index_type extent
[GFC_MAX_DIMENSIONS
];
255 index_type stride
[GFC_MAX_DIMENSIONS
];
264 if (x
->dim
[0].stride
== 0)
265 x
->dim
[0].stride
= 1;
267 dim
= GFC_DESCRIPTOR_RANK (x
);
269 for (n
= 0; n
< dim
; n
++)
272 stride
[n
] = x
->dim
[n
].stride
;
273 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
280 __gthread_mutex_lock (&random_lock
);
284 /* random_r8 (dest); */
285 kiss
= ((GFC_UINTEGER_8
)kiss_random_kernel ()) << 32;
286 kiss
+= kiss_random_kernel ();
287 *dest
= normalize_r8_i8 (kiss
, ~(GFC_UINTEGER_8
) 0);
289 /* Advance to the next element. */
292 /* Advance to the next source element. */
294 while (count
[n
] == extent
[n
])
296 /* When we get to the end of a dimension, reset it and increment
297 the next dimension. */
299 /* We could precalculate these products, but this is a less
300 frequently used path so probably not worth it. */
301 dest
-= stride
[n
] * extent
[n
];
315 __gthread_mutex_unlock (&random_lock
);
318 /* random_seed is used to seed the PRNG with either a default
319 set of seeds or user specified set of seeds. random_seed
320 must be called with no argument or exactly one argument. */
323 random_seed (GFC_INTEGER_4
*size
, gfc_array_i4
*put
, gfc_array_i4
*get
)
327 __gthread_mutex_lock (&random_lock
);
329 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
331 /* From the standard: "If no argument is present, the processor assigns
332 a processor-dependent value to the seed." */
333 kiss_seed
[0] = kiss_default_seed
[0];
334 kiss_seed
[1] = kiss_default_seed
[1];
335 kiss_seed
[2] = kiss_default_seed
[2];
336 kiss_seed
[3] = kiss_default_seed
[3];
344 /* If the rank of the array is not 1, abort. */
345 if (GFC_DESCRIPTOR_RANK (put
) != 1)
346 runtime_error ("Array rank of PUT is not 1.");
348 /* If the array is too small, abort. */
349 if (((put
->dim
[0].ubound
+ 1 - put
->dim
[0].lbound
)) < kiss_size
)
350 runtime_error ("Array size of PUT is too small.");
352 if (put
->dim
[0].stride
== 0)
353 put
->dim
[0].stride
= 1;
355 /* This code now should do correct strides. */
356 for (i
= 0; i
< kiss_size
; i
++)
357 kiss_seed
[i
] =(GFC_UINTEGER_4
) put
->data
[i
* put
->dim
[0].stride
];
360 /* Return the seed to GET data. */
363 /* If the rank of the array is not 1, abort. */
364 if (GFC_DESCRIPTOR_RANK (get
) != 1)
365 runtime_error ("Array rank of GET is not 1.");
367 /* If the array is too small, abort. */
368 if (((get
->dim
[0].ubound
+ 1 - get
->dim
[0].lbound
)) < kiss_size
)
369 runtime_error ("Array size of GET is too small.");
371 if (get
->dim
[0].stride
== 0)
372 get
->dim
[0].stride
= 1;
374 /* This code now should do correct strides. */
375 for (i
= 0; i
< kiss_size
; i
++)
376 get
->data
[i
* get
->dim
[0].stride
] = (GFC_INTEGER_4
) kiss_seed
[i
];
379 __gthread_mutex_unlock (&random_lock
);
381 iexport(random_seed
);
384 #ifndef __GTHREAD_MUTEX_INIT
385 static void __attribute__((constructor
))
388 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock
);