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., 59 Temple Place - Suite 330,
30 Boston, MA 02111-1307, USA. */
32 #include "libgfortran.h"
34 extern void random_r4 (GFC_REAL_4 *);
35 iexport_proto(random_r4);
37 extern void random_r8 (GFC_REAL_8 *);
38 iexport_proto(random_r8);
40 extern void arandom_r4 (gfc_array_r4 *);
41 export_proto(arandom_r4);
43 extern void arandom_r8 (gfc_array_r8 *);
44 export_proto(arandom_r8);
48 /* The Mersenne Twister code is currently commented out due to
50 (1) Simple user specified seeds lead to really bad sequences for
51 nearly 100000 random numbers.
52 (2) open(), read(), and close() are not properly declared via header
54 (3) The global index i is abused and causes unexpected behavior with
58 The algorithm was taken from the paper :
60 Mersenne Twister: 623-dimensionally equidistributed
61 uniform pseudorandom generator.
66 Which appeared in the: ACM Transactions on Modelling and Computer
67 Simulations: Special Issue on Uniform Random Number
68 Generation. ( Early in 1998 ). */
73 #include <sys/types.h>
81 /*Use the 'big' generator by default ( period -> 2**19937 ). */
85 /* Define the necessary constants for the algorithm. */
90 N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
92 #define M_A 0x9908B0DF
93 #define T_B 0x9D2C5680
94 #define T_C 0xEFC60000
98 N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17
100 #define M_A 0xE4BD75F5
101 #define T_B 0x655E5280
102 #define T_C 0xFFD58000
106 static unsigned int seed[N];
108 /* This is the routine which handles the seeding of the generator,
109 and also reading and writing of the seed. */
112 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
114 /* Initialize the seed in system dependent manner. */
115 if (get == NULL && put == NULL && size == NULL)
118 fd = open ("/dev/urandom", O_RDONLY);
121 /* We dont have urandom. */
122 GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed;
123 for (i = 0; i < N; i++)
125 s = s * 29943829 - 1;
131 /* Using urandom, might have a length issue. */
132 read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N);
138 /* Return the size of the seed */
145 /* if we have gotten to this pount we have a get or put
146 * now we check it the array fulfills the demands in the standard .
149 /* Set the seed to PUT data */
152 /* if the rank of the array is not 1 abort */
153 if (GFC_DESCRIPTOR_RANK (put) != 1)
156 /* if the array is too small abort */
157 if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N)
160 /* If this is the case the array is a temporary */
161 if (put->dim[0].stride == 0)
164 /* This code now should do correct strides. */
165 for (i = 0; i < N; i++)
166 seed[i] = put->data[i * put->dim[0].stride];
169 /* Return the seed to GET data */
172 /* if the rank of the array is not 1 abort */
173 if (GFC_DESCRIPTOR_RANK (get) != 1)
176 /* if the array is too small abort */
177 if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N)
180 /* If this is the case the array is a temporary */
181 if (get->dim[0].stride == 0)
184 /* This code now should do correct strides. */
185 for (i = 0; i < N; i++)
186 get->data[i * get->dim[0].stride] = seed[i];
189 iexport(random_seed);
191 /* Here is the internal routine which generates the random numbers
192 in 'batches' based upon the need for a new batch.
193 It's an integer based routine known as 'Mersenne Twister'.
194 This implementation still lacks 'tempering' and a good verification,
195 but gives very good metrics. */
198 random_generate (void)
203 /* Generate batch of N. */
205 for (k = 0, m = M; k < N - 1; k++)
207 y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1));
208 seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
213 y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1));
214 seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A);
218 /* A routine to return a REAL(KIND=4). */
221 random_r4 (GFC_REAL_4 * harv)
223 /* Regenerate if we need to. */
227 /* Convert uint32 to REAL(KIND=4). */
228 *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
229 (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
233 /* A routine to return a REAL(KIND=8). */
236 random_r8 (GFC_REAL_8 * harv)
238 /* Regenerate if we need to, may waste one 32-bit value. */
242 /* Convert two uint32 to a REAL(KIND=8). */
243 *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
244 (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
249 /* Code to handle arrays will follow here. */
251 /* REAL(KIND=4) REAL array. */
254 arandom_r4 (gfc_array_r4 * harv)
256 index_type count[GFC_MAX_DIMENSIONS];
257 index_type extent[GFC_MAX_DIMENSIONS];
258 index_type stride[GFC_MAX_DIMENSIONS];
266 if (harv->dim[0].stride == 0)
267 harv->dim[0].stride = 1;
269 dim = GFC_DESCRIPTOR_RANK (harv);
271 for (n = 0; n < dim; n++)
274 stride[n] = harv->dim[n].stride;
275 extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
284 /* Set the elements. */
286 /* regenerate if we need to */
290 /* Convert uint32 to float in a hopefully g95 compiant manner */
291 *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] /
292 (GFC_REAL_4) (~(GFC_UINTEGER_4) 0));
294 /* Advance to the next element. */
297 /* Advance to the next source element. */
299 while (count[n] == extent[n])
301 /* When we get to the end of a dimension,
302 reset it and increment
303 the next dimension. */
305 /* We could precalculate these products,
307 frequently used path so proabably not worth it. */
308 dest -= stride[n] * extent[n];
324 /* REAL(KIND=8) array. */
327 arandom_r8 (gfc_array_r8 * harv)
329 index_type count[GFC_MAX_DIMENSIONS];
330 index_type extent[GFC_MAX_DIMENSIONS];
331 index_type stride[GFC_MAX_DIMENSIONS];
339 if (harv->dim[0].stride == 0)
340 harv->dim[0].stride = 1;
342 dim = GFC_DESCRIPTOR_RANK (harv);
344 for (n = 0; n < dim; n++)
347 stride[n] = harv->dim[n].stride;
348 extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound;
357 /* Set the elements. */
359 /* regenerate if we need to, may waste one 32-bit value */
363 /* Convert two uint32 to a REAL(KIND=8). */
364 *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) /
365 (GFC_REAL_8) (~(GFC_UINTEGER_8) 0);
368 /* Advance to the next element. */
371 /* Advance to the next source element. */
373 while (count[n] == extent[n])
375 /* When we get to the end of a dimension,
376 reset it and increment
377 the next dimension. */
379 /* We could precalculate these products,
381 frequently used path so proabably not worth it. */
382 dest -= stride[n] * extent[n];
400 /* George Marsaglia's KISS (Keep It Simple Stupid) random number generator.
404 (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
406 (2) A 3-shift shift-register generator with a period of 2^32-1,
407 (3) Two 16-bit multiply-with-carry generators with a period of
408 597273182964842497 > 2^59.
410 The overall period exceeds 2^123.
412 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
414 The above web site has an archive of a newsgroup posting from George
415 Marsaglia with the statement:
417 Subject: Random numbers for C: Improvements.
418 Date: Fri, 15 Jan 1999 11:41:47 -0500
419 From: George Marsaglia <geo@stat.fsu.edu>
420 Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
421 References: <369B5E30.65A55FD1@stat.fsu.edu>
422 Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
425 As I hoped, several suggestions have led to
426 improvements in the code for RNG's I proposed for
427 use in C. (See the thread "Random numbers for C: Some
428 suggestions" in previous postings.) The improved code
431 A question of copyright has also been raised. Unlike
432 DIEHARD, there is no copyright on the code below. You
433 are free to use it in any way you want, but you may
434 wish to acknowledge the source, as a courtesy.
436 "There is no copyright on the code below." included the original
439 #define GFC_SL(k, n) ((k)^((k)<<(n)))
440 #define GFC_SR(k, n) ((k)^((k)>>(n)))
442 static const GFC_INTEGER_4 kiss_size = 4;
443 #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069}
444 static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED;
445 static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED;
447 /* kiss_random_kernel() returns an integer value in the range of
448 (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
449 should be uniform. */
451 static GFC_UINTEGER_4
452 kiss_random_kernel(void)
456 kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885;
457 kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5);
458 kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16);
459 kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16);
460 kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3];
465 /* This function produces a REAL(4) value from the uniform distribution
469 random_r4 (GFC_REAL_4 *x)
473 kiss = kiss_random_kernel ();
474 /* Burn a random number, so the REAL*4 and REAL*8 functions
475 produce similar sequences of random numbers. */
476 kiss_random_kernel ();
477 *x = normalize_r4_i4 (kiss, ~(GFC_UINTEGER_4) 0);
481 /* This function produces a REAL(8) value from the uniform distribution
485 random_r8 (GFC_REAL_8 *x)
489 kiss = ((GFC_UINTEGER_8)kiss_random_kernel ()) << 32;
490 kiss += kiss_random_kernel ();
491 *x = normalize_r8_i8 (kiss, ~(GFC_UINTEGER_8) 0);
495 /* This function fills a REAL(4) array with values from the uniform
496 distribution with range [0,1). */
499 arandom_r4 (gfc_array_r4 *x)
501 index_type count[GFC_MAX_DIMENSIONS];
502 index_type extent[GFC_MAX_DIMENSIONS];
503 index_type stride[GFC_MAX_DIMENSIONS];
511 if (x->dim[0].stride == 0)
512 x->dim[0].stride = 1;
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;
531 /* Advance to the next element. */
534 /* Advance to the next source element. */
536 while (count[n] == extent[n])
538 /* When we get to the end of a dimension, reset it and increment
539 the next dimension. */
541 /* We could precalculate these products, but this is a less
542 frequently used path so probably not worth it. */
543 dest -= stride[n] * extent[n];
559 /* This function fills a REAL(8) array with values from the uniform
560 distribution with range [0,1). */
563 arandom_r8 (gfc_array_r8 *x)
565 index_type count[GFC_MAX_DIMENSIONS];
566 index_type extent[GFC_MAX_DIMENSIONS];
567 index_type stride[GFC_MAX_DIMENSIONS];
575 if (x->dim[0].stride == 0)
576 x->dim[0].stride = 1;
578 dim = GFC_DESCRIPTOR_RANK (x);
580 for (n = 0; n < dim; n++)
583 stride[n] = x->dim[n].stride;
584 extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound;
595 /* Advance to the next element. */
598 /* Advance to the next source element. */
600 while (count[n] == extent[n])
602 /* When we get to the end of a dimension, reset it and increment
603 the next dimension. */
605 /* We could precalculate these products, but this is a less
606 frequently used path so probably not worth it. */
607 dest -= stride[n] * extent[n];
623 /* random_seed is used to seed the PRNG with either a default
624 set of seeds or user specified set of seeds. random_seed
625 must be called with no argument or exactly one argument. */
628 random_seed (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get)
632 if (size == NULL && put == NULL && get == NULL)
634 /* From the standard: "If no argument is present, the processor assigns
635 a processor-dependent value to the seed." */
636 kiss_seed[0] = kiss_default_seed[0];
637 kiss_seed[1] = kiss_default_seed[1];
638 kiss_seed[2] = kiss_default_seed[2];
639 kiss_seed[3] = kiss_default_seed[3];
647 /* If the rank of the array is not 1, abort. */
648 if (GFC_DESCRIPTOR_RANK (put) != 1)
649 runtime_error ("Array rank of PUT is not 1.");
651 /* If the array is too small, abort. */
652 if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
653 runtime_error ("Array size of PUT is too small.");
655 if (put->dim[0].stride == 0)
656 put->dim[0].stride = 1;
658 /* This code now should do correct strides. */
659 for (i = 0; i < kiss_size; i++)
660 kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
663 /* Return the seed to GET data. */
666 /* If the rank of the array is not 1, abort. */
667 if (GFC_DESCRIPTOR_RANK (get) != 1)
668 runtime_error ("Array rank of GET is not 1.");
670 /* If the array is too small, abort. */
671 if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
672 runtime_error ("Array size of GET is too small.");
674 if (get->dim[0].stride == 0)
675 get->dim[0].stride = 1;
677 /* This code now should do correct strides. */
678 for (i = 0; i < kiss_size; i++)
679 get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
682 iexport(random_seed);
684 #endif /* mersenne twister */