X-Git-Url: http://git.sourceforge.jp/view?p=pf3gnuchains%2Fgcc-fork.git;a=blobdiff_plain;f=libgfortran%2Fintrinsics%2Frandom.c;h=f64f21c5b80d12956a226aa0f2160eeaabc7c5cd;hp=73a6cede0f0cd550545a7ae83d4b38decb1a2f59;hb=3d3b790db32448aed311a1df5c182f639ee7ec97;hpb=b0fb613e02eb095d8c4b9655404c29e46f8bf495 diff --git a/libgfortran/intrinsics/random.c b/libgfortran/intrinsics/random.c index 73a6cede0f0..f64f21c5b80 100644 --- a/libgfortran/intrinsics/random.c +++ b/libgfortran/intrinsics/random.c @@ -1,276 +1,399 @@ /* Implementation of the RANDOM intrinsics - Copyright 2002, 2004 Free Software Foundation, Inc. + Copyright 2002, 2004, 2005, 2006, 2007 Free Software Foundation, Inc. Contributed by Lars Segerlund and Steve Kargl. This file is part of the GNU Fortran 95 runtime library (libgfortran). Libgfortran is free software; you can redistribute it and/or -modify it under the terms of the GNU Lesser General Public +modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either -version 2.1 of the License, or (at your option) any later version. +version 2 of the License, or (at your option) any later version. + +In addition to the permissions in the GNU General Public License, the +Free Software Foundation gives you unlimited permission to link the +compiled version of this file into combinations with other programs, +and to distribute those combinations without any restriction coming +from the use of this file. (The General Public License restrictions +do apply in other respects; for example, they cover modification of +the file, and distribution when not linked into a combine +executable.) Ligbfortran is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -GNU Lesser General Public License for more details. +GNU General Public License for more details. -You should have received a copy of the GNU Lesser General Public -License along with libgfor; see the file COPYING.LIB. If not, -write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, -Boston, MA 02111-1307, USA. */ +You should have received a copy of the GNU General Public +License along with libgfortran; see the file COPYING. If not, +write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. */ -#if 0 +#include "config.h" +#include "libgfortran.h" +#include +#include -/* The Mersenne Twister code is currently commented out due to +extern void random_r4 (GFC_REAL_4 *); +iexport_proto(random_r4); - (1) Simple user specified seeds lead to really bad sequences for - nearly 100000 random numbers. - (2) open(), read(), and close() are not properly declared via header - files. - (3) The global index i is abused and causes unexpected behavior with - GET and PUT. - (4) See PR 15619. +extern void random_r8 (GFC_REAL_8 *); +iexport_proto(random_r8); - The algorithm was taken from the paper : +extern void arandom_r4 (gfc_array_r4 *); +export_proto(arandom_r4); - Mersenne Twister: 623-dimensionally equidistributed - uniform pseudorandom generator. +extern void arandom_r8 (gfc_array_r8 *); +export_proto(arandom_r8); - by: Makoto Matsumoto - Takuji Nishimura +#ifdef HAVE_GFC_REAL_10 - Which appeared in the: ACM Transactions on Modelling and Computer - Simulations: Special Issue on Uniform Random Number - Generation. ( Early in 1998 ). */ +extern void random_r10 (GFC_REAL_10 *); +iexport_proto(random_r10); +extern void arandom_r10 (gfc_array_r10 *); +export_proto(arandom_r10); -#include "config.h" -#include -#include -#include -#include -#include - -#ifdef HAVE_UNISTD_H -#include #endif -#include "libgfortran.h" +#ifdef HAVE_GFC_REAL_16 + +extern void random_r16 (GFC_REAL_16 *); +iexport_proto(random_r16); -/*Use the 'big' generator by default ( period -> 2**19937 ). */ +extern void arandom_r16 (gfc_array_r16 *); +export_proto(arandom_r16); -#define MT19937 +#endif -/* Define the necessary constants for the algorithm. */ +#ifdef __GTHREAD_MUTEX_INIT +static __gthread_mutex_t random_lock = __GTHREAD_MUTEX_INIT; +#else +static __gthread_mutex_t random_lock; +#endif -#ifdef MT19937 -enum constants +/* Helper routines to map a GFC_UINTEGER_* to the corresponding + GFC_REAL_* types in the range of [0,1). If GFC_REAL_*_RADIX are 2 + or 16, respectively, we mask off the bits that don't fit into the + correct GFC_REAL_*, convert to the real type, then multiply by the + correct offset. +*/ + + +static inline void +rnumber_4 (GFC_REAL_4 *f, GFC_UINTEGER_4 v) { - N = 624, M = 397, R = 19, TU = 11, TS = 7, TT = 15, TL = 17 -}; -#define M_A 0x9908B0DF -#define T_B 0x9D2C5680 -#define T_C 0xEFC60000 + GFC_UINTEGER_4 mask; +#if GFC_REAL_4_RADIX == 2 + mask = ~ (GFC_UINTEGER_4) 0u << (32 - GFC_REAL_4_DIGITS); +#elif GFC_REAL_4_RADIX == 16 + mask = ~ (GFC_UINTEGER_4) 0u << ((8 - GFC_REAL_4_DIGITS) * 4); #else -enum constants +#error "GFC_REAL_4_RADIX has unknown value" +#endif + v = v & mask; + *f = (GFC_REAL_4) v * (GFC_REAL_4) 0x1.p-32f; +} + +static inline void +rnumber_8 (GFC_REAL_8 *f, GFC_UINTEGER_8 v) { - N = 351, M = 175, R = 19, TU = 11, TS = 7, TT = 15, TL = 17 -}; -#define M_A 0xE4BD75F5 -#define T_B 0x655E5280 -#define T_C 0xFFD58000 + GFC_UINTEGER_8 mask; +#if GFC_REAL_8_RADIX == 2 + mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_8_DIGITS); +#elif GFC_REAL_8_RADIX == 16 + mask = ~ (GFC_UINTEGER_8) 0u << (16 - GFC_REAL_8_DIGITS) * 4); +#else +#error "GFC_REAL_8_RADIX has unknown value" #endif + v = v & mask; + *f = (GFC_REAL_8) v * (GFC_REAL_8) 0x1.p-64; +} -static int i = N; -static unsigned int seed[N]; +#ifdef HAVE_GFC_REAL_10 -/* This is the routine which handles the seeding of the generator, - and also reading and writing of the seed. */ +static inline void +rnumber_10 (GFC_REAL_10 *f, GFC_UINTEGER_8 v) +{ + GFC_UINTEGER_8 mask; +#if GFC_REAL_10_RADIX == 2 + mask = ~ (GFC_UINTEGER_8) 0u << (64 - GFC_REAL_10_DIGITS); +#elif GFC_REAL_10_RADIX == 16 + mask = ~ (GFC_UINTEGER_10) 0u << ((16 - GFC_REAL_10_DIGITS) * 4); +#else +#error "GFC_REAL_10_RADIX has unknown value" +#endif + v = v & mask; + *f = (GFC_REAL_10) v * (GFC_REAL_10) 0x1.p-64; +} +#endif -void -random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put, - const gfc_array_i4 * get) +#ifdef HAVE_GFC_REAL_16 + +/* For REAL(KIND=16), we only need to mask off the lower bits. */ + +static inline void +rnumber_16 (GFC_REAL_16 *f, GFC_UINTEGER_8 v1, GFC_UINTEGER_8 v2) { - /* Initialize the seed in system dependent manner. */ - if (get == NULL && put == NULL && size == NULL) - { - int fd; - fd = open ("/dev/urandom", O_RDONLY); - if (fd == 0) - { - /* We dont have urandom. */ - GFC_UINTEGER_4 s = (GFC_UINTEGER_4) seed; - for (i = 0; i < N; i++) - { - s = s * 29943829 - 1; - seed[i] = s; - } - } - else - { - /* Using urandom, might have a length issue. */ - read (fd, &seed[0], sizeof (GFC_UINTEGER_4) * N); - close (fd); - } - return; - } + GFC_UINTEGER_8 mask; +#if GFC_REAL_16_RADIX == 2 + mask = ~ (GFC_UINTEGER_8) 0u << (128 - GFC_REAL_16_DIGITS); +#elif GFC_REAL_16_RADIX == 16 + mask = ~ (GFC_UINTEGER_8) 0u << ((32 - GFC_REAL_16_DIGITS) * 4); +#else +#error "GFC_REAL_16_RADIX has unknown value" +#endif + v2 = v2 & mask; + *f = (GFC_REAL_16) v1 * (GFC_REAL_16) 0x1.p-64 + + (GFC_REAL_16) v2 * (GFC_REAL_16) 0x1.p-128; +} +#endif +/* libgfortran previously had a Mersenne Twister, taken from the paper: + + Mersenne Twister: 623-dimensionally equidistributed + uniform pseudorandom generator. - /* Return the size of the seed */ - if (size != NULL) - { - *size = N; - return; - } + by Makoto Matsumoto & Takuji Nishimura + which appeared in the: ACM Transactions on Modelling and Computer + Simulations: Special Issue on Uniform Random Number + Generation. ( Early in 1998 ). - /* if we have gotten to this pount we have a get or put - * now we check it the array fulfills the demands in the standard . - */ + The Mersenne Twister code was replaced due to - /* Set the seed to PUT data */ - if (put != NULL) - { - /* if the rank of the array is not 1 abort */ - if (GFC_DESCRIPTOR_RANK (put) != 1) - abort (); + (1) Simple user specified seeds lead to really bad sequences for + nearly 100000 random numbers. + (2) open(), read(), and close() were not properly declared via header + files. + (3) The global index i was abused and caused unexpected behavior with + GET and PUT. + (4) See PR 15619. - /* if the array is too small abort */ - if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < N) - abort (); - /* If this is the case the array is a temporary */ - if (put->dim[0].stride == 0) - return; + libgfortran currently uses George Marsaglia's KISS (Keep It Simple Stupid) + random number generator. This PRNG combines: - /* This code now should do correct strides. */ - for (i = 0; i < N; i++) - seed[i] = put->data[i * put->dim[0].stride]; - } + (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period + of 2^32, + (2) A 3-shift shift-register generator with a period of 2^32-1, + (3) Two 16-bit multiply-with-carry generators with a period of + 597273182964842497 > 2^59. - /* Return the seed to GET data */ - if (get != NULL) - { - /* if the rank of the array is not 1 abort */ - if (GFC_DESCRIPTOR_RANK (get) != 1) - abort (); + The overall period exceeds 2^123. - /* if the array is too small abort */ - if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < N) - abort (); + http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu - /* If this is the case the array is a temporary */ - if (get->dim[0].stride == 0) - return; + The above web site has an archive of a newsgroup posting from George + Marsaglia with the statement: - /* This code now should do correct strides. */ - for (i = 0; i < N; i++) - get->data[i * get->dim[0].stride] = seed[i]; - } + Subject: Random numbers for C: Improvements. + Date: Fri, 15 Jan 1999 11:41:47 -0500 + From: George Marsaglia + Message-ID: <369F6FCA.74C7C041@stat.fsu.edu> + References: <369B5E30.65A55FD1@stat.fsu.edu> + Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis + Lines: 93 + + As I hoped, several suggestions have led to + improvements in the code for RNG's I proposed for + use in C. (See the thread "Random numbers for C: Some + suggestions" in previous postings.) The improved code + is listed below. + + A question of copyright has also been raised. Unlike + DIEHARD, there is no copyright on the code below. You + are free to use it in any way you want, but you may + wish to acknowledge the source, as a courtesy. + +"There is no copyright on the code below." included the original +KISS algorithm. */ + +/* We use three KISS random number generators, with different + seeds. + As a matter of Quality of Implementation, the random numbers + we generate for different REAL kinds, starting from the same + seed, are always the same up to the precision of these types. + We do this by using three generators with different seeds, the + first one always for the most significant bits, the second one + for bits 33..64 (if present in the REAL kind), and the third one + (called twice) for REAL(16). +*/ + +#define GFC_SL(k, n) ((k)^((k)<<(n))) +#define GFC_SR(k, n) ((k)^((k)>>(n))) + +/* Reference for the seed: + From: "George Marsaglia" + Newsgroups: sci.math + Message-ID: + + The KISS RNG uses four seeds, x, y, z, c, + with 0<=x<2^32, 0> 16); + seed[3] = 30903 * (seed[3] & 65535) + (seed[3] >> 16); + kiss = seed[0] + seed[1] + (seed[2] << 16) + seed[3]; + + return kiss; } -/* Here is the internal routine which generates the random numbers - in 'batches' based upon the need for a new batch. - It's an integer based routine known as 'Mersenne Twister'. - This implementation still lacks 'tempering' and a good verification, - but gives very good metrics. */ +/* This function produces a REAL(4) value from the uniform distribution + with range [0,1). */ -static void -random_generate (void) +void +random_r4 (GFC_REAL_4 *x) { - /* 32 bits. */ - GFC_UINTEGER_4 y; + GFC_UINTEGER_4 kiss; - /* Generate batch of N. */ - int k, m; - for (k = 0, m = M; k < N - 1; k++) - { - y = (seed[k] & (-1 << R)) | (seed[k + 1] & ((1u << R) - 1)); - seed[k] = seed[m] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A); - if (++m >= N) - m = 0; - } + __gthread_mutex_lock (&random_lock); + kiss = kiss_random_kernel (kiss_seed_1); + rnumber_4 (x, kiss); + __gthread_mutex_unlock (&random_lock); +} +iexport(random_r4); + +/* This function produces a REAL(8) value from the uniform distribution + with range [0,1). */ + +void +random_r8 (GFC_REAL_8 *x) +{ + GFC_UINTEGER_8 kiss; - y = (seed[N - 1] & (-1 << R)) | (seed[0] & ((1u << R) - 1)); - seed[N - 1] = seed[M - 1] ^ (y >> 1) ^ (-(GFC_INTEGER_4) (y & 1) & M_A); - i = 0; + __gthread_mutex_lock (&random_lock); + kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; + kiss += kiss_random_kernel (kiss_seed_2); + rnumber_8 (x, kiss); + __gthread_mutex_unlock (&random_lock); } +iexport(random_r8); -/* A routine to return a REAL(KIND=4). */ +#ifdef HAVE_GFC_REAL_10 + +/* This function produces a REAL(10) value from the uniform distribution + with range [0,1). */ -#define random_r4 prefix(random_r4) void -random_r4 (GFC_REAL_4 * harv) +random_r10 (GFC_REAL_10 *x) { - /* Regenerate if we need to. */ - if (i >= N) - random_generate (); + GFC_UINTEGER_8 kiss; - /* Convert uint32 to REAL(KIND=4). */ - *harv = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] / - (GFC_REAL_4) (~(GFC_UINTEGER_4) 0)); + __gthread_mutex_lock (&random_lock); + kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; + kiss += kiss_random_kernel (kiss_seed_2); + rnumber_10 (x, kiss); + __gthread_mutex_unlock (&random_lock); } +iexport(random_r10); + +#endif -/* A routine to return a REAL(KIND=8). */ +/* This function produces a REAL(16) value from the uniform distribution + with range [0,1). */ + +#ifdef HAVE_GFC_REAL_16 -#define random_r8 prefix(random_r8) void -random_r8 (GFC_REAL_8 * harv) +random_r16 (GFC_REAL_16 *x) { - /* Regenerate if we need to, may waste one 32-bit value. */ - if ((i + 1) >= N) - random_generate (); - - /* Convert two uint32 to a REAL(KIND=8). */ - *harv = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) / - (GFC_REAL_8) (~(GFC_UINTEGER_8) 0); - i += 2; + GFC_UINTEGER_8 kiss1, kiss2; + + __gthread_mutex_lock (&random_lock); + kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; + kiss1 += kiss_random_kernel (kiss_seed_2); + + kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32; + kiss2 += kiss_random_kernel (kiss_seed_3); + + rnumber_16 (x, kiss1, kiss2); + __gthread_mutex_unlock (&random_lock); } +iexport(random_r16); -/* Code to handle arrays will follow here. */ -/* REAL(KIND=4) REAL array. */ +#endif +/* This function fills a REAL(4) array with values from the uniform + distribution with range [0,1). */ -#define arandom_r4 prefix(arandom_r4) void -arandom_r4 (gfc_array_r4 * harv) +arandom_r4 (gfc_array_r4 *x) { - index_type count[GFC_MAX_DIMENSIONS - 1]; - index_type extent[GFC_MAX_DIMENSIONS - 1]; - index_type stride[GFC_MAX_DIMENSIONS - 1]; + index_type count[GFC_MAX_DIMENSIONS]; + index_type extent[GFC_MAX_DIMENSIONS]; + index_type stride[GFC_MAX_DIMENSIONS]; index_type stride0; index_type dim; GFC_REAL_4 *dest; + GFC_UINTEGER_4 kiss; int n; - dest = harv->data; - - if (harv->dim[0].stride == 0) - harv->dim[0].stride = 1; + dest = x->data; - dim = GFC_DESCRIPTOR_RANK (harv); + dim = GFC_DESCRIPTOR_RANK (x); for (n = 0; n < dim; n++) { count[n] = 0; - stride[n] = harv->dim[n].stride; - extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound; + stride[n] = x->dim[n].stride; + extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound; if (extent[n] <= 0) - return; + return; } stride0 = stride[0]; + __gthread_mutex_lock (&random_lock); + while (dest) { - /* Set the elements. */ - - /* regenerate if we need to */ - if (i >= N) - random_generate (); - - /* Convert uint32 to float in a hopefully g95 compiant manner */ - *dest = (GFC_REAL_4) ((GFC_REAL_4) (GFC_UINTEGER_4) seed[i++] / - (GFC_REAL_4) (~(GFC_UINTEGER_4) 0)); + /* random_r4 (dest); */ + kiss = kiss_random_kernel (kiss_seed_1); + rnumber_4 (dest, kiss); /* Advance to the next element. */ dest += stride0; @@ -278,74 +401,67 @@ arandom_r4 (gfc_array_r4 * harv) /* Advance to the next source element. */ n = 0; while (count[n] == extent[n]) - { - /* When we get to the end of a dimension, - reset it and increment - the next dimension. */ - count[n] = 0; - /* We could precalculate these products, - but this is a less - frequently used path so proabably not worth it. */ - dest -= stride[n] * extent[n]; - n++; - if (n == dim) - { - dest = NULL; - break; - } - else - { - count[n]++; - dest += stride[n]; - } - } + { + /* When we get to the end of a dimension, reset it and increment + the next dimension. */ + count[n] = 0; + /* We could precalculate these products, but this is a less + frequently used path so probably not worth it. */ + dest -= stride[n] * extent[n]; + n++; + if (n == dim) + { + dest = NULL; + break; + } + else + { + count[n]++; + dest += stride[n]; + } + } } + __gthread_mutex_unlock (&random_lock); } -/* REAL(KIND=8) array. */ +/* This function fills a REAL(8) array with values from the uniform + distribution with range [0,1). */ -#define arandom_r8 prefix(arandom_r8) void -arandom_r8 (gfc_array_r8 * harv) +arandom_r8 (gfc_array_r8 *x) { - index_type count[GFC_MAX_DIMENSIONS - 1]; - index_type extent[GFC_MAX_DIMENSIONS - 1]; - index_type stride[GFC_MAX_DIMENSIONS - 1]; + index_type count[GFC_MAX_DIMENSIONS]; + index_type extent[GFC_MAX_DIMENSIONS]; + index_type stride[GFC_MAX_DIMENSIONS]; index_type stride0; index_type dim; GFC_REAL_8 *dest; + GFC_UINTEGER_8 kiss; int n; - dest = harv->data; - - if (harv->dim[0].stride == 0) - harv->dim[0].stride = 1; + dest = x->data; - dim = GFC_DESCRIPTOR_RANK (harv); + dim = GFC_DESCRIPTOR_RANK (x); for (n = 0; n < dim; n++) { count[n] = 0; - stride[n] = harv->dim[n].stride; - extent[n] = harv->dim[n].ubound + 1 - harv->dim[n].lbound; + stride[n] = x->dim[n].stride; + extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound; if (extent[n] <= 0) - return; + return; } stride0 = stride[0]; + __gthread_mutex_lock (&random_lock); + while (dest) { - /* Set the elements. */ - - /* regenerate if we need to, may waste one 32-bit value */ - if ((i + 1) >= N) - random_generate (); - - /* Convert two uint32 to a REAL(KIND=8). */ - *dest = ((GFC_REAL_8) ((((GFC_UINTEGER_8) seed[i+1]) << 32) + seed[i])) / - (GFC_REAL_8) (~(GFC_UINTEGER_8) 0); - i += 2; + /* random_r8 (dest); */ + kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; + kiss += kiss_random_kernel (kiss_seed_2); + rnumber_8 (dest, kiss); /* Advance to the next element. */ dest += stride0; @@ -353,162 +469,48 @@ arandom_r8 (gfc_array_r8 * harv) /* Advance to the next source element. */ n = 0; while (count[n] == extent[n]) - { - /* When we get to the end of a dimension, - reset it and increment - the next dimension. */ - count[n] = 0; - /* We could precalculate these products, - but this is a less - frequently used path so proabably not worth it. */ - dest -= stride[n] * extent[n]; - n++; - if (n == dim) - { - dest = NULL; - break; - } - else - { - count[n]++; - dest += stride[n]; - } - } - } -} -#endif /* Mersenne Twister code */ - - -/* George Marsaglia's KISS (Keep It Simple Stupid) random number generator. - - This PRNG combines: - - (1) The congruential generator x(n)=69069*x(n-1)+1327217885 with a period - of 2^32, - (2) A 3-shift shift-register generator with a period of 2^32-1, - (3) Two 16-bit multiply-with-carry generators with a period of - 597273182964842497 > 2^59. - - The overall period exceeds 2^123. - - http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369F6FCA.74C7C041@stat.fsu.edu - - The above web site has an archive of a newsgroup posting from George - Marsaglia with the statement: - - Subject: Random numbers for C: Improvements. - Date: Fri, 15 Jan 1999 11:41:47 -0500 - From: George Marsaglia - Message-ID: <369F6FCA.74C7C041@stat.fsu.edu> - References: <369B5E30.65A55FD1@stat.fsu.edu> - Newsgroups: sci.stat.math,sci.math,sci.math.numer-analysis - Lines: 93 - - As I hoped, several suggestions have led to - improvements in the code for RNG's I proposed for - use in C. (See the thread "Random numbers for C: Some - suggestions" in previous postings.) The improved code - is listed below. - - A question of copyright has also been raised. Unlike - DIEHARD, there is no copyright on the code below. You - are free to use it in any way you want, but you may - wish to acknowledge the source, as a courtesy. - -"There is no copyright on the code below." included the original -KISS algorithm. */ - -#include "config.h" -#include "libgfortran.h" - -#define GFC_SL(k, n) ((k)^((k)<<(n))) -#define GFC_SR(k, n) ((k)^((k)>>(n))) - -static const GFC_INTEGER_4 kiss_size = 4; -#define KISS_DEFAULT_SEED {123456789, 362436069, 521288629, 916191069}; -static const GFC_UINTEGER_4 kiss_default_seed[4] = KISS_DEFAULT_SEED; -static GFC_UINTEGER_4 kiss_seed[4] = KISS_DEFAULT_SEED; - -/* kiss_random_kernel() returns an integer value in the range of - (0, GFC_UINTEGER_4_HUGE]. The distribution of pseudorandom numbers - should be uniform. */ - -static GFC_UINTEGER_4 -kiss_random_kernel(void) -{ - - GFC_UINTEGER_4 kiss; - - kiss_seed[0] = 69069 * kiss_seed[0] + 1327217885; - kiss_seed[1] = GFC_SL(GFC_SR(GFC_SL(kiss_seed[1],13),17),5); - kiss_seed[2] = 18000 * (kiss_seed[2] & 65535) + (kiss_seed[2] >> 16); - kiss_seed[3] = 30903 * (kiss_seed[3] & 65535) + (kiss_seed[3] >> 16); - kiss = kiss_seed[0] + kiss_seed[1] + (kiss_seed[2] << 16) + kiss_seed[3]; - - return kiss; - -} - -/* This function produces a REAL(4) value in the uniform distribution - with range [0,1). */ - -void -prefix(random_r4) (GFC_REAL_4 *x) -{ - - GFC_UINTEGER_4 kiss; - - do - { - kiss = kiss_random_kernel (); - *x = (GFC_REAL_4)kiss / (GFC_REAL_4)(~(GFC_UINTEGER_4) 0); - /* Burn a random number, so the REAL*4 and REAL*8 functions - produce similar sequences of random numbers. */ - kiss = kiss_random_kernel (); + { + /* When we get to the end of a dimension, reset it and increment + the next dimension. */ + count[n] = 0; + /* We could precalculate these products, but this is a less + frequently used path so probably not worth it. */ + dest -= stride[n] * extent[n]; + n++; + if (n == dim) + { + dest = NULL; + break; + } + else + { + count[n]++; + dest += stride[n]; + } + } } - while (*x == 1.0); - + __gthread_mutex_unlock (&random_lock); } -/* This function produces a REAL(8) value from the uniform distribution - with range [0,1). */ - -void -prefix(random_r8) (GFC_REAL_8 *x) -{ - - GFC_UINTEGER_8 kiss; - - do - { - kiss = (((GFC_UINTEGER_8)kiss_random_kernel ()) << 32) - + kiss_random_kernel (); - *x = (GFC_REAL_8)kiss / (GFC_REAL_8)(~(GFC_UINTEGER_8) 0); - } - while (*x == 0); +#ifdef HAVE_GFC_REAL_10 -} - -/* This function fills a REAL(4) array with values from the uniform +/* This function fills a REAL(10) array with values from the uniform distribution with range [0,1). */ void -prefix(arandom_r4) (gfc_array_r4 *x) +arandom_r10 (gfc_array_r10 *x) { - - index_type count[GFC_MAX_DIMENSIONS - 1]; - index_type extent[GFC_MAX_DIMENSIONS - 1]; - index_type stride[GFC_MAX_DIMENSIONS - 1]; + index_type count[GFC_MAX_DIMENSIONS]; + index_type extent[GFC_MAX_DIMENSIONS]; + index_type stride[GFC_MAX_DIMENSIONS]; index_type stride0; index_type dim; - GFC_REAL_4 *dest; + GFC_REAL_10 *dest; + GFC_UINTEGER_8 kiss; int n; dest = x->data; - if (x->dim[0].stride == 0) - x->dim[0].stride = 1; - dim = GFC_DESCRIPTOR_RANK (x); for (n = 0; n < dim; n++) @@ -522,9 +524,14 @@ prefix(arandom_r4) (gfc_array_r4 *x) stride0 = stride[0]; + __gthread_mutex_lock (&random_lock); + while (dest) { - prefix(random_r4) (dest); + /* random_r10 (dest); */ + kiss = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; + kiss += kiss_random_kernel (kiss_seed_2); + rnumber_10 (dest, kiss); /* Advance to the next element. */ dest += stride0; @@ -552,28 +559,30 @@ prefix(arandom_r4) (gfc_array_r4 *x) } } } + __gthread_mutex_unlock (&random_lock); } -/* This function fills a REAL(8) array with valuse from the uniform +#endif + +#ifdef HAVE_GFC_REAL_16 + +/* This function fills a REAL(16) array with values from the uniform distribution with range [0,1). */ void -prefix(arandom_r8) (gfc_array_r8 *x) +arandom_r16 (gfc_array_r16 *x) { - - index_type count[GFC_MAX_DIMENSIONS - 1]; - index_type extent[GFC_MAX_DIMENSIONS - 1]; - index_type stride[GFC_MAX_DIMENSIONS - 1]; + index_type count[GFC_MAX_DIMENSIONS]; + index_type extent[GFC_MAX_DIMENSIONS]; + index_type stride[GFC_MAX_DIMENSIONS]; index_type stride0; index_type dim; - GFC_REAL_8 *dest; + GFC_REAL_16 *dest; + GFC_UINTEGER_8 kiss1, kiss2; int n; dest = x->data; - if (x->dim[0].stride == 0) - x->dim[0].stride = 1; - dim = GFC_DESCRIPTOR_RANK (x); for (n = 0; n < dim; n++) @@ -587,9 +596,18 @@ prefix(arandom_r8) (gfc_array_r8 *x) stride0 = stride[0]; + __gthread_mutex_lock (&random_lock); + while (dest) { - prefix(random_r8) (dest); + /* random_r16 (dest); */ + kiss1 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_1)) << 32; + kiss1 += kiss_random_kernel (kiss_seed_2); + + kiss2 = ((GFC_UINTEGER_8) kiss_random_kernel (kiss_seed_3)) << 32; + kiss2 += kiss_random_kernel (kiss_seed_3); + + rnumber_16 (dest, kiss1, kiss2); /* Advance to the next element. */ dest += stride0; @@ -617,68 +635,133 @@ prefix(arandom_r8) (gfc_array_r8 *x) } } } + __gthread_mutex_unlock (&random_lock); } -/* prefix(random_seed) is used to seed the PRNG with either a default - set of seeds or user specified set of seed. prefix(random_seed) +#endif + +/* random_seed is used to seed the PRNG with either a default + set of seeds or user specified set of seeds. random_seed must be called with no argument or exactly one argument. */ void -random_seed (GFC_INTEGER_4 *size, gfc_array_i4 * put, - gfc_array_i4 * get) +random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) { - int i; + __gthread_mutex_lock (&random_lock); + + /* Check that we only have one argument present. */ + if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) + runtime_error ("RANDOM_SEED should have at most one argument present."); + + /* From the standard: "If no argument is present, the processor assigns + a processor-dependent value to the seed." */ if (size == NULL && put == NULL && get == NULL) - { - /* From the standard: "If no argument is present, the processor assigns - a processor-dependent value to the seed." */ - kiss_seed[0] = kiss_default_seed[0]; - kiss_seed[1] = kiss_default_seed[1]; - kiss_seed[2] = kiss_default_seed[2]; - kiss_seed[3] = kiss_default_seed[3]; - } + for (i = 0; i < kiss_size; i++) + kiss_seed[i] = kiss_default_seed[i]; if (size != NULL) *size = kiss_size; if (put != NULL) { - /* If the rank of the array is not 1, abort */ + /* If the rank of the array is not 1, abort. */ if (GFC_DESCRIPTOR_RANK (put) != 1) runtime_error ("Array rank of PUT is not 1."); - /* If the array is too small, abort */ + /* If the array is too small, abort. */ if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size) runtime_error ("Array size of PUT is too small."); - if (put->dim[0].stride == 0) - put->dim[0].stride = 1; - - /* This code now should do correct strides. */ + /* This code now should do correct strides. */ for (i = 0; i < kiss_size; i++) - kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride]; + kiss_seed[i] = (GFC_UINTEGER_4) put->data[i * put->dim[0].stride]; } - /* Return the seed to GET data */ + /* Return the seed to GET data. */ if (get != NULL) { - /* If the rank of the array is not 1, abort. */ + /* If the rank of the array is not 1, abort. */ if (GFC_DESCRIPTOR_RANK (get) != 1) runtime_error ("Array rank of GET is not 1."); - /* If the array is too small, abort. */ + /* If the array is too small, abort. */ if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size) runtime_error ("Array size of GET is too small."); - if (get->dim[0].stride == 0) - get->dim[0].stride = 1; - - /* This code now should do correct strides. */ + /* This code now should do correct strides. */ for (i = 0; i < kiss_size; i++) get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i]; } + + __gthread_mutex_unlock (&random_lock); +} +iexport(random_seed_i4); + + +void +random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get) +{ + int i; + + __gthread_mutex_lock (&random_lock); + + /* Check that we only have one argument present. */ + if ((size ? 1 : 0) + (put ? 1 : 0) + (get ? 1 : 0) > 1) + runtime_error ("RANDOM_SEED should have at most one argument present."); + + /* From the standard: "If no argument is present, the processor assigns + a processor-dependent value to the seed." */ + if (size == NULL && put == NULL && get == NULL) + for (i = 0; i < kiss_size; i++) + kiss_seed[i] = kiss_default_seed[i]; + + if (size != NULL) + *size = kiss_size / 2; + + if (put != NULL) + { + /* If the rank of the array is not 1, abort. */ + if (GFC_DESCRIPTOR_RANK (put) != 1) + runtime_error ("Array rank of PUT is not 1."); + + /* If the array is too small, abort. */ + if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size / 2) + runtime_error ("Array size of PUT is too small."); + + /* This code now should do correct strides. */ + for (i = 0; i < kiss_size; i += 2) + memcpy (&kiss_seed[i], &(put->data[i * put->dim[0].stride]), + sizeof (GFC_UINTEGER_8)); + } + + /* Return the seed to GET data. */ + if (get != NULL) + { + /* If the rank of the array is not 1, abort. */ + if (GFC_DESCRIPTOR_RANK (get) != 1) + runtime_error ("Array rank of GET is not 1."); + + /* If the array is too small, abort. */ + if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size / 2) + runtime_error ("Array size of GET is too small."); + + /* This code now should do correct strides. */ + for (i = 0; i < kiss_size; i += 2) + memcpy (&(get->data[i * get->dim[0].stride]), &kiss_seed[i], + sizeof (GFC_UINTEGER_8)); + } + + __gthread_mutex_unlock (&random_lock); } +iexport(random_seed_i8); +#ifndef __GTHREAD_MUTEX_INIT +static void __attribute__((constructor)) +init (void) +{ + __GTHREAD_MUTEX_INIT_FUNCTION (&random_lock); +} +#endif