From f36ac248b3556cd7f13c76bc19a0300afcdd3ea8 Mon Sep 17 00:00:00 2001 From: pbrook Date: Sun, 30 May 2004 10:49:50 +0000 Subject: [PATCH] * iresolve.c (gfc_resolve_random_number): Clean up conditional. libgfortran/ * libgfortran.h (random_seed): Update prototype. * intrinsics/random.c: Disable old implementation and add new one. testsuite/ * gfortran.fortran-torture/execute/random_1.f90: New test. git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@82443 138bc75d-0d04-0410-961f-82ee72b054a4 --- gcc/fortran/ChangeLog | 4 + gcc/fortran/iresolve.c | 9 +- gcc/testsuite/ChangeLog | 4 + .../gfortran.fortran-torture/execute/random_1.f90 | 33 ++ libgfortran/ChangeLog | 5 + libgfortran/intrinsics/random.c | 340 ++++++++++++++++++++- libgfortran/libgfortran.h | 4 +- 7 files changed, 382 insertions(+), 17 deletions(-) create mode 100644 gcc/testsuite/gfortran.fortran-torture/execute/random_1.f90 diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog index ac795569089..3bc18098771 100644 --- a/gcc/fortran/ChangeLog +++ b/gcc/fortran/ChangeLog @@ -1,3 +1,7 @@ +2004-05-30 Steven G. Kargl + + * iresolve.c (gfc_resolve_random_number): Clean up conditional. + 2004-05-29 Steven G. Kargl * simplify.c (gfc_simplify_log): Remove useless line of code. diff --git a/gcc/fortran/iresolve.c b/gcc/fortran/iresolve.c index 46e38037f60..f1da7329c48 100644 --- a/gcc/fortran/iresolve.c +++ b/gcc/fortran/iresolve.c @@ -1363,12 +1363,15 @@ gfc_resolve_random_number (gfc_code * c ATTRIBUTE_UNUSED) int kind; kind = c->ext.actual->expr->ts.kind; - name = gfc_get_string ((c->ext.actual->expr->rank == 0) ? - PREFIX("random_r%d") : PREFIX("arandom_r%d"), - kind); + if (c->ext.actual->expr->rank == 0) + name = gfc_get_string (PREFIX("random_r%d"), kind); + else + name = gfc_get_string (PREFIX("arandom_r%d"), kind); + c->resolved_sym = gfc_get_intrinsic_sub_symbol (name); } + /* Determine if the arguments to SYSTEM_CLOCK are INTEGER(4) or INTEGER(8) */ void diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog index ed4955f7fe5..0535cfe5ce3 100644 --- a/gcc/testsuite/ChangeLog +++ b/gcc/testsuite/ChangeLog @@ -1,3 +1,7 @@ +2004-05-30 Steven G. Kargl + + * gfortran.fortran-torture/execute/random_1.f90: New test. + 2004-05-28 Ziemowit Laski * g++.dg/ext/altivec-10.C: New test. diff --git a/gcc/testsuite/gfortran.fortran-torture/execute/random_1.f90 b/gcc/testsuite/gfortran.fortran-torture/execute/random_1.f90 new file mode 100644 index 00000000000..900a724da6c --- /dev/null +++ b/gcc/testsuite/gfortran.fortran-torture/execute/random_1.f90 @@ -0,0 +1,33 @@ +! PR15619 +! Check that random_seed works as expected. +! Does not check the quality of random numbers, hence should never fail. +program test_random + implicit none + integer, allocatable :: seed(:) + real, dimension(10) :: a, b + integer n; + + call random_seed (size=n) + allocate (seed(n)) + + ! Exercise the generator a bit. + call random_number (a) + + ! Remeber the seed and get 10 more. + call random_seed (get=seed) + call random_number (a) + + ! Get the same 10 numbers in two blocks, remebering the seed in the middle + call random_seed (put=seed) + call random_number (b(1:5)) + call random_seed(get=seed) + call random_number (b(6:10)) + if (any (a .ne. b)) call abort + + ! Get the last 5 numbers again. + call random_seed (put=seed) + call random_number (b(6:10)) + if (any (a .ne. b)) call abort +end program + + diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog index 92fc985e211..def98c62670 100644 --- a/libgfortran/ChangeLog +++ b/libgfortran/ChangeLog @@ -1,3 +1,8 @@ +2004-05-30 Steven G. Kargl + + * libgfortran.h (random_seed): Update prototype. + * intrinsics/random.c: Disable old implementation and add new one. + 2004-05-30 Andreas Jaeger * intrinsics/random.c: Include unistd.h for close and read diff --git a/libgfortran/intrinsics/random.c b/libgfortran/intrinsics/random.c index c1539243c0e..bfda3437f91 100644 --- a/libgfortran/intrinsics/random.c +++ b/libgfortran/intrinsics/random.c @@ -1,18 +1,7 @@ /* Implementation of the RANDOM intrinsics Copyright 2002, 2004 Free Software Foundation, Inc. Contributed by Lars Segerlund - - The algorithm was taken from the paper : - - Mersenne Twister: 623-dimensionally equidistributed - uniform pseudorandom generator. - - 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 ). + and Steve Kargl. This file is part of the GNU Fortran 95 runtime library (libgfortran). @@ -31,6 +20,31 @@ 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. */ +#if 0 + +/* The Mersenne Twister code is currently commented out due to + + (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. + + The algorithm was taken from the paper : + + Mersenne Twister: 623-dimensionally equidistributed + uniform pseudorandom generator. + + 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 ). */ + + #include "config.h" #include #include @@ -362,4 +376,306 @@ arandom_r8 (gfc_array_r8 * harv) } } } +#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); + } + while (*x == 1.0); + +} + +/* 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); + +} + +/* This function fills a REAL(4) array with values from the uniform + distribution with range [0,1). */ + +void +prefix(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 stride0; + index_type dim; + GFC_REAL_4 *dest; + 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++) + { + count[n] = 0; + stride[n] = x->dim[n].stride; + extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound; + if (extent[n] <= 0) + return; + } + + stride0 = stride[0]; + + while (dest) + { + prefix(random_r4) (dest); + + /* Advance to the next element. */ + dest += stride0; + count[0]++; + /* 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 probably not worth it. */ + dest -= stride[n] * extent[n]; + n++; + if (n == dim) + { + dest = NULL; + break; + } + else + { + count[n]++; + dest += stride[n]; + } + } + } +} + +/* This function fills a REAL(8) array with valuse from the uniform + distribution with range [0,1). */ + +void +prefix(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 stride0; + index_type dim; + GFC_REAL_8 *dest; + 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++) + { + count[n] = 0; + stride[n] = x->dim[n].stride; + extent[n] = x->dim[n].ubound + 1 - x->dim[n].lbound; + if (extent[n] <= 0) + return; + } + + stride0 = stride[0]; + + while (dest) + { + prefix(random_r8) (dest); + + /* Advance to the next element. */ + dest += stride0; + count[0]++; + /* 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 probably not worth it. */ + dest -= stride[n] * extent[n]; + n++; + if (n == dim) + { + dest = NULL; + break; + } + else + { + count[n]++; + dest += stride[n]; + } + } + } +} + +/* 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) + 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) +{ + + int i; + + 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]; + } + + if (size != NULL) + *size = kiss_size; + + 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) + 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. */ + for (i = 0; i < kiss_size; i++) + kiss_seed[i] =(GFC_UINTEGER_4) put->data[i * put->dim[0].stride]; + } + + /* 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) + 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. */ + for (i = 0; i < kiss_size; i++) + get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i]; + } +} + diff --git a/libgfortran/libgfortran.h b/libgfortran/libgfortran.h index 6234e25cb72..a4c27597280 100644 --- a/libgfortran/libgfortran.h +++ b/libgfortran/libgfortran.h @@ -399,8 +399,8 @@ GFC_INTEGER_4 compare_string (GFC_INTEGER_4, const char *, /* random.c */ #define random_seed prefix(random_seed) -void random_seed (GFC_INTEGER_4 * size, const gfc_array_i4 * put, - const gfc_array_i4 * get); +void random_seed (GFC_INTEGER_4 * size, gfc_array_i4 * put, + gfc_array_i4 * get); #endif -- 2.11.0