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
;
56 /* The Mersenne Twister code is currently commented out due to
58 (1) Simple user specified seeds lead to really bad sequences for
59 nearly 100000 random numbers.
60 (2) open(), read(), and close() are not properly declared via header
62 (3) The global index i is abused and causes unexpected behavior with
66 The algorithm was taken from the paper :
68 Mersenne Twister: 623-dimensionally equidistributed
69 uniform pseudorandom generator.
74 Which appeared in the: ACM Transactions on Modelling and Computer
75 Simulations: Special Issue on Uniform Random Number
76 Generation. ( Early in 1998 ). */
81 #include <sys/types.h>
89 /*Use the 'big' generator by default ( period -> 2**19937 ). */
93 /* Define the necessary constants for the algorithm. */
98 N
= 624, M
= 397, R
= 19, TU
= 11, TS
= 7, TT
= 15, TL
= 17
100 #define M_A 0x9908B0DF
101 #define T_B 0x9D2C5680
102 #define T_C 0xEFC60000
106 N
= 351, M
= 175, R
= 19, TU
= 11, TS
= 7, TT
= 15, TL
= 17
108 #define M_A 0xE4BD75F5
109 #define T_B 0x655E5280
110 #define T_C 0xFFD58000
114 static unsigned int seed
[N
];
116 /* This is the routine which handles the seeding of the generator,
117 and also reading and writing of the seed. */
120 random_seed (GFC_INTEGER_4
*size
, gfc_array_i4
*put
, gfc_array_i4
*get
)
122 __gthread_mutex_lock (&random_lock
);
124 /* Initialize the seed in system dependent manner. */
125 if (get
== NULL
&& put
== NULL
&& size
== NULL
)
128 fd
= open ("/dev/urandom", O_RDONLY
);
131 /* We dont have urandom. */
132 GFC_UINTEGER_4 s
= (GFC_UINTEGER_4
) seed
;
133 for (i
= 0; i
< N
; i
++)
135 s
= s
* 29943829 - 1;
141 /* Using urandom, might have a length issue. */
142 read (fd
, &seed
[0], sizeof (GFC_UINTEGER_4
) * N
);
149 /* Return the size of the seed */
156 /* if we have gotten to this pount we have a get or put
157 * now we check it the array fulfills the demands in the standard .
160 /* Set the seed to PUT data */
163 /* if the rank of the array is not 1 abort */
164 if (GFC_DESCRIPTOR_RANK (put
) != 1)
167 /* if the array is too small abort */
168 if (((put
->dim
[0].ubound
+ 1 - put
->dim
[0].lbound
)) < N
)
171 /* If this is the case the array is a temporary */
172 if (put
->dim
[0].stride
== 0)
175 /* This code now should do correct strides. */
176 for (i
= 0; i
< N
; i
++)
177 seed
[i
] = put
->data
[i
* put
->dim
[0].stride
];
180 /* Return the seed to GET data */
183 /* if the rank of the array is not 1 abort */
184 if (GFC_DESCRIPTOR_RANK (get
) != 1)
187 /* if the array is too small abort */
188 if (((get
->dim
[0].ubound
+ 1 - get
->dim
[0].lbound
)) < N
)
191 /* If this is the case the array is a temporary */
192 if (get
->dim
[0].stride
== 0)
195 /* This code now should do correct strides. */
196 for (i
= 0; i
< N
; i
++)
197 get
->data
[i
* get
->dim
[0].stride
] = seed
[i
];
201 __gthread_mutex_unlock (&random_lock
);
203 iexport(random_seed
);
205 /* Here is the internal routine which generates the random numbers
206 in 'batches' based upon the need for a new batch.
207 It's an integer based routine known as 'Mersenne Twister'.
208 This implementation still lacks 'tempering' and a good verification,
209 but gives very good metrics. */
212 random_generate (void)
217 /* Generate batch of N. */
219 for (k
= 0, m
= M
; k
< N
- 1; k
++)
221 y
= (seed
[k
] & (-1 << R
)) | (seed
[k
+ 1] & ((1u << R
) - 1));
222 seed
[k
] = seed
[m
] ^ (y
>> 1) ^ (-(GFC_INTEGER_4
) (y
& 1) & M_A
);
227 y
= (seed
[N
- 1] & (-1 << R
)) | (seed
[0] & ((1u << R
) - 1));
228 seed
[N
- 1] = seed
[M
- 1] ^ (y
>> 1) ^ (-(GFC_INTEGER_4
) (y
& 1) & M_A
);
232 /* A routine to return a REAL(KIND=4). */
235 random_r4 (GFC_REAL_4
* harv
)
237 __gthread_mutex_lock (&random_lock
);
239 /* Regenerate if we need to. */
243 /* Convert uint32 to REAL(KIND=4). */
244 *harv
= (GFC_REAL_4
) ((GFC_REAL_4
) (GFC_UINTEGER_4
) seed
[i
++] /
245 (GFC_REAL_4
) (~(GFC_UINTEGER_4
) 0));
246 __gthread_mutex_unlock (&random_lock
);
250 /* A routine to return a REAL(KIND=8). */
253 random_r8 (GFC_REAL_8
* harv
)
255 __gthread_mutex_lock (&random_lock
);
257 /* Regenerate if we need to, may waste one 32-bit value. */
261 /* Convert two uint32 to a REAL(KIND=8). */
262 *harv
= ((GFC_REAL_8
) ((((GFC_UINTEGER_8
) seed
[i
+1]) << 32) + seed
[i
])) /
263 (GFC_REAL_8
) (~(GFC_UINTEGER_8
) 0);
265 __gthread_mutex_unlock (&random_lock
);
269 /* Code to handle arrays will follow here. */
271 /* REAL(KIND=4) REAL array. */
274 arandom_r4 (gfc_array_r4
* harv
)
276 index_type count
[GFC_MAX_DIMENSIONS
];
277 index_type extent
[GFC_MAX_DIMENSIONS
];
278 index_type stride
[GFC_MAX_DIMENSIONS
];
286 if (harv
->dim
[0].stride
== 0)
287 harv
->dim
[0].stride
= 1;
289 dim
= GFC_DESCRIPTOR_RANK (harv
);
291 for (n
= 0; n
< dim
; n
++)
294 stride
[n
] = harv
->dim
[n
].stride
;
295 extent
[n
] = harv
->dim
[n
].ubound
+ 1 - harv
->dim
[n
].lbound
;
302 __gthread_mutex_lock (&random_lock
);
306 /* Set the elements. */
308 /* regenerate if we need to */
312 /* Convert uint32 to float in a hopefully g95 compiant manner */
313 *dest
= (GFC_REAL_4
) ((GFC_REAL_4
) (GFC_UINTEGER_4
) seed
[i
++] /
314 (GFC_REAL_4
) (~(GFC_UINTEGER_4
) 0));
316 /* Advance to the next element. */
319 /* Advance to the next source element. */
321 while (count
[n
] == extent
[n
])
323 /* When we get to the end of a dimension,
324 reset it and increment
325 the next dimension. */
327 /* We could precalculate these products,
329 frequently used path so proabably not worth it. */
330 dest
-= stride
[n
] * extent
[n
];
345 __gthread_mutex_unlock (&random_lock
);
348 /* REAL(KIND=8) array. */
351 arandom_r8 (gfc_array_r8
* harv
)
353 index_type count
[GFC_MAX_DIMENSIONS
];
354 index_type extent
[GFC_MAX_DIMENSIONS
];
355 index_type stride
[GFC_MAX_DIMENSIONS
];
363 if (harv
->dim
[0].stride
== 0)
364 harv
->dim
[0].stride
= 1;
366 dim
= GFC_DESCRIPTOR_RANK (harv
);
368 for (n
= 0; n
< dim
; n
++)
371 stride
[n
] = harv
->dim
[n
].stride
;
372 extent
[n
] = harv
->dim
[n
].ubound
+ 1 - harv
->dim
[n
].lbound
;
379 __gthread_mutex_lock (&random_lock
);
383 /* Set the elements. */
385 /* regenerate if we need to, may waste one 32-bit value */
389 /* Convert two uint32 to a REAL(KIND=8). */
390 *dest
= ((GFC_REAL_8
) ((((GFC_UINTEGER_8
) seed
[i
+1]) << 32) + seed
[i
])) /
391 (GFC_REAL_8
) (~(GFC_UINTEGER_8
) 0);
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,
402 reset it and increment
403 the next dimension. */
405 /* We could precalculate these products,
407 frequently used path so proabably not worth it. */
408 dest
-= stride
[n
] * extent
[n
];
423 __gthread_mutex_unlock (&random_lock
);
428 /* George Marsaglia's KISS (Keep It Simple Stupid) random number generator.
432 (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period
434 (2) A 3-shift shift-register generator with a period of 2^32-1,
435 (3) Two 16-bit multiply-with-carry generators with a period of
436 597273182964842497 > 2^59.
438 The overall period exceeds 2^123.
440 http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu
442 The above web site has an archive of a newsgroup posting from George
443 Marsaglia with the statement:
445 Subject: Random numbers for C: Improvements.
446 Date: Fri, 15 Jan 1999 11:41:47 -0500
447 From: George Marsaglia <geo@stat.fsu.edu>
448 Message-ID: <369F6FCA.74C7C041@stat.fsu.edu>
449 References: <369B5E30.65A55FD1@stat.fsu.edu>
450 Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis
453 As I hoped, several suggestions have led to
454 improvements in the code for RNG's I proposed for
455 use in C. (See the thread "Random numbers for C: Some
456 suggestions" in previous postings.) The improved code
459 A question of copyright has also been raised. Unlike
460 DIEHARD, there is no copyright on the code below. You
461 are free to use it in any way you want, but you may
462 wish to acknowledge the source, as a courtesy.
464 "There is no copyright on the code below." included the original
467 #define GFC_SL(k, n) ((k)^((k)<<(n)))
468 #define GFC_SR(k, n) ((k)^((k)>>(n)))
470 static const GFC_INTEGER_4 kiss_size
= 4;
471 #define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069}
472 static const GFC_UINTEGER_4 kiss_default_seed
[4] = KISS_DEFAULT_SEED
;
473 static GFC_UINTEGER_4 kiss_seed
[4] = KISS_DEFAULT_SEED
;
475 /* kiss_random_kernel() returns an integer value in the range of
476 (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers
477 should be uniform. */
479 static GFC_UINTEGER_4
480 kiss_random_kernel(void)
484 kiss_seed
[0] = 69069 * kiss_seed
[0] + 1327217885;
485 kiss_seed
[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed
[1],13),17),5);
486 kiss_seed
[2] = 18000 * (kiss_seed
[2] & 65535) + (kiss_seed
[2] >> 16);
487 kiss_seed
[3] = 30903 * (kiss_seed
[3] & 65535) + (kiss_seed
[3] >> 16);
488 kiss
= kiss_seed
[0] + kiss_seed
[1] + (kiss_seed
[2] << 16) + kiss_seed
[3];
493 /* This function produces a REAL(4) value from the uniform distribution
497 random_r4 (GFC_REAL_4
*x
)
501 __gthread_mutex_lock (&random_lock
);
502 kiss
= kiss_random_kernel ();
503 /* Burn a random number, so the REAL*4 and REAL*8 functions
504 produce similar sequences of random numbers. */
505 kiss_random_kernel ();
506 *x
= normalize_r4_i4 (kiss
, ~(GFC_UINTEGER_4
) 0);
507 __gthread_mutex_unlock (&random_lock
);
511 /* This function produces a REAL(8) value from the uniform distribution
515 random_r8 (GFC_REAL_8
*x
)
519 __gthread_mutex_lock (&random_lock
);
520 kiss
= ((GFC_UINTEGER_8
)kiss_random_kernel ()) << 32;
521 kiss
+= kiss_random_kernel ();
522 *x
= normalize_r8_i8 (kiss
, ~(GFC_UINTEGER_8
) 0);
523 __gthread_mutex_unlock (&random_lock
);
527 /* This function fills a REAL(4) array with values from the uniform
528 distribution with range [0,1). */
531 arandom_r4 (gfc_array_r4
*x
)
533 index_type count
[GFC_MAX_DIMENSIONS
];
534 index_type extent
[GFC_MAX_DIMENSIONS
];
535 index_type stride
[GFC_MAX_DIMENSIONS
];
544 if (x
->dim
[0].stride
== 0)
545 x
->dim
[0].stride
= 1;
547 dim
= GFC_DESCRIPTOR_RANK (x
);
549 for (n
= 0; n
< dim
; n
++)
552 stride
[n
] = x
->dim
[n
].stride
;
553 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
560 __gthread_mutex_lock (&random_lock
);
564 /* random_r4 (dest); */
565 kiss
= kiss_random_kernel ();
566 /* Burn a random number, so the REAL*4 and REAL*8 functions
567 produce similar sequences of random numbers. */
568 kiss_random_kernel ();
569 *dest
= normalize_r4_i4 (kiss
, ~(GFC_UINTEGER_4
) 0);
571 /* Advance to the next element. */
574 /* Advance to the next source element. */
576 while (count
[n
] == extent
[n
])
578 /* When we get to the end of a dimension, reset it and increment
579 the next dimension. */
581 /* We could precalculate these products, but this is a less
582 frequently used path so probably not worth it. */
583 dest
-= stride
[n
] * extent
[n
];
597 __gthread_mutex_unlock (&random_lock
);
600 /* This function fills a REAL(8) array with values from the uniform
601 distribution with range [0,1). */
604 arandom_r8 (gfc_array_r8
*x
)
606 index_type count
[GFC_MAX_DIMENSIONS
];
607 index_type extent
[GFC_MAX_DIMENSIONS
];
608 index_type stride
[GFC_MAX_DIMENSIONS
];
617 if (x
->dim
[0].stride
== 0)
618 x
->dim
[0].stride
= 1;
620 dim
= GFC_DESCRIPTOR_RANK (x
);
622 for (n
= 0; n
< dim
; n
++)
625 stride
[n
] = x
->dim
[n
].stride
;
626 extent
[n
] = x
->dim
[n
].ubound
+ 1 - x
->dim
[n
].lbound
;
633 __gthread_mutex_lock (&random_lock
);
637 /* random_r8 (dest); */
638 kiss
= ((GFC_UINTEGER_8
)kiss_random_kernel ()) << 32;
639 kiss
+= kiss_random_kernel ();
640 *dest
= normalize_r8_i8 (kiss
, ~(GFC_UINTEGER_8
) 0);
642 /* Advance to the next element. */
645 /* Advance to the next source element. */
647 while (count
[n
] == extent
[n
])
649 /* When we get to the end of a dimension, reset it and increment
650 the next dimension. */
652 /* We could precalculate these products, but this is a less
653 frequently used path so probably not worth it. */
654 dest
-= stride
[n
] * extent
[n
];
668 __gthread_mutex_unlock (&random_lock
);
671 /* random_seed is used to seed the PRNG with either a default
672 set of seeds or user specified set of seeds. random_seed
673 must be called with no argument or exactly one argument. */
676 random_seed (GFC_INTEGER_4
*size
, gfc_array_i4
*put
, gfc_array_i4
*get
)
680 __gthread_mutex_lock (&random_lock
);
682 if (size
== NULL
&& put
== NULL
&& get
== NULL
)
684 /* From the standard: "If no argument is present, the processor assigns
685 a processor-dependent value to the seed." */
686 kiss_seed
[0] = kiss_default_seed
[0];
687 kiss_seed
[1] = kiss_default_seed
[1];
688 kiss_seed
[2] = kiss_default_seed
[2];
689 kiss_seed
[3] = kiss_default_seed
[3];
697 /* If the rank of the array is not 1, abort. */
698 if (GFC_DESCRIPTOR_RANK (put
) != 1)
699 runtime_error ("Array rank of PUT is not 1.");
701 /* If the array is too small, abort. */
702 if (((put
->dim
[0].ubound
+ 1 - put
->dim
[0].lbound
)) < kiss_size
)
703 runtime_error ("Array size of PUT is too small.");
705 if (put
->dim
[0].stride
== 0)
706 put
->dim
[0].stride
= 1;
708 /* This code now should do correct strides. */
709 for (i
= 0; i
< kiss_size
; i
++)
710 kiss_seed
[i
] =(GFC_UINTEGER_4
) put
->data
[i
* put
->dim
[0].stride
];
713 /* Return the seed to GET data. */
716 /* If the rank of the array is not 1, abort. */
717 if (GFC_DESCRIPTOR_RANK (get
) != 1)
718 runtime_error ("Array rank of GET is not 1.");
720 /* If the array is too small, abort. */
721 if (((get
->dim
[0].ubound
+ 1 - get
->dim
[0].lbound
)) < kiss_size
)
722 runtime_error ("Array size of GET is too small.");
724 if (get
->dim
[0].stride
== 0)
725 get
->dim
[0].stride
= 1;
727 /* This code now should do correct strides. */
728 for (i
= 0; i
< kiss_size
; i
++)
729 get
->data
[i
* get
->dim
[0].stride
] = (GFC_INTEGER_4
) kiss_seed
[i
];
732 __gthread_mutex_unlock (&random_lock
);
734 iexport(random_seed
);
736 #endif /* mersenne twister */
738 #ifndef __GTHREAD_MUTEX_INIT
739 static void __attribute__((constructor
))
742 __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock
);