All, Here is a potentially contentious patch for PR fortran/52879.
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=52879 As demonstrated in that PR, it is possible to provide RANDOM_SEED with sets of seeds that result in a very poor sequences of PRN. Gfortran's RANDOM_NUMBER uses 3 independent KISS PRNG to generate the bits of the significands of the real PRN. Each KISS PRNG requires 4 seeds. This patch removes the need for a user to specify 12 seeds. It uses a Lehmer linear congruent generator to generate the 12 KISS seeds. This LCG requires a single seed. I have selected to have seed(1)=0 correspond to the current default set of KISS seeds. Internally, the LCG declares the seed as a GFC_UINTEGER_8 (ie., uint64_t) type, so for gfortran's default integer kind one has 2**31-1 possible seeds and if one use -fdefault-integer-8 then one has 2**32 possible seeds. It is possible to extend the patch to allow a skip count such that for example seed(2)=10 would use every tenth value in the LCG sequence. I haven't pursued this as it appears to be an unneeded complication. The only testsuite failure is for gfortran.dg/random_seed_1.f90, which is a test designed to specifically test for the 12 KISS seeds. This test will be removed if the patch is OK. -- Steve Index: gcc/fortran/check.c =================================================================== --- gcc/fortran/check.c (revision 207633) +++ gcc/fortran/check.c (working copy) @@ -4925,16 +4925,9 @@ gfc_check_random_number (gfc_expr *harve bool gfc_check_random_seed (gfc_expr *size, gfc_expr *put, gfc_expr *get) { - unsigned int nargs = 0, kiss_size; + unsigned int nargs = 0; locus *where = NULL; mpz_t put_size, get_size; - bool have_gfc_real_16; /* Try and mimic HAVE_GFC_REAL_16 in libgfortran. */ - - have_gfc_real_16 = gfc_validate_kind (BT_REAL, 16, true) != -1; - - /* Keep the number of bytes in sync with kiss_size in - libgfortran/intrinsics/random.c. */ - kiss_size = (have_gfc_real_16 ? 48 : 32) / gfc_default_integer_kind; if (size != NULL) { @@ -4975,13 +4968,6 @@ gfc_check_random_seed (gfc_expr *size, g if (!kind_value_check (put, 1, gfc_default_integer_kind)) return false; - - if (gfc_array_size (put, &put_size) - && mpz_get_ui (put_size) < kiss_size) - gfc_error ("Size of '%s' argument of '%s' intrinsic at %L " - "too small (%i/%i)", - gfc_current_intrinsic_arg[1]->name, gfc_current_intrinsic, - where, (int) mpz_get_ui (put_size), kiss_size); } if (get != NULL) @@ -5007,13 +4993,6 @@ gfc_check_random_seed (gfc_expr *size, g if (!kind_value_check (get, 2, gfc_default_integer_kind)) return false; - - if (gfc_array_size (get, &get_size) - && mpz_get_ui (get_size) < kiss_size) - gfc_error ("Size of '%s' argument of '%s' intrinsic at %L " - "too small (%i/%i)", - gfc_current_intrinsic_arg[2]->name, gfc_current_intrinsic, - where, (int) mpz_get_ui (get_size), kiss_size); } /* RANDOM_SEED may not have more than one non-optional argument. */ Index: libgfortran/intrinsics/random.c =================================================================== --- libgfortran/intrinsics/random.c (revision 207633) +++ libgfortran/intrinsics/random.c (working copy) @@ -224,9 +224,20 @@ KISS algorithm. */ z=0,c=0 and z=2^32-1,c=698769068 should be avoided. */ -/* Any modifications to the seeds that change kiss_size below need to be - reflected in check.c (gfc_check_random_seed) to enable correct - compile-time checking of PUT size for the RANDOM_SEED intrinsic. */ +/* The 3 KISS generators require 3 sets for 4 seeds. It is possible for a + a user to specify a set of seeds that is inadequate for reseeding the + generators. For example, change only the value of one of the 12 seeds + may not reset the PRNG sequences. To work-around this issue, use a + simple Lehmer linear congruential generator that requires only a single + seed to generate the 12 KISS seeds. */ + +static const GFC_INTEGER_4 lehmer_lcg_size = 1; +static GFC_UINTEGER_8 lehmer_lcg_seed = 0; +static GFC_UINTEGER_4 +lehmer_lcg(GFC_UINTEGER_4 a) +{ + return ((GFC_UINTEGER_8)a * 279470273UL) % 4294967291UL; +} #define KISS_DEFAULT_SEED_1 123456789, 362436069, 521288629, 316191069 #define KISS_DEFAULT_SEED_2 987654321, 458629013, 582859209, 438195021 @@ -636,28 +647,6 @@ arandom_r16 (gfc_array_r16 *x) #endif - -static void -scramble_seed (unsigned char *dest, unsigned char *src, int size) -{ - int i; - - for (i = 0; i < size; i++) - dest[(i % 2) * (size / 2) + i / 2] = src[i]; -} - - -static void -unscramble_seed (unsigned char *dest, unsigned char *src, int size) -{ - int i; - - for (i = 0; i < size; i++) - dest[i] = src[(i % 2) * (size / 2) + i / 2]; -} - - - /* 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. */ @@ -666,7 +655,7 @@ void random_seed_i4 (GFC_INTEGER_4 *size, gfc_array_i4 *put, gfc_array_i4 *get) { int i; - unsigned char seed[4*kiss_size]; + GFC_UINTEGER_4 seed; __gthread_mutex_lock (&random_lock); @@ -676,54 +665,33 @@ random_seed_i4 (GFC_INTEGER_4 *size, gfc /* 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) + if (size == NULL && put == NULL && get == NULL) { + lehmer_lcg_seed = 0; for (i = 0; i < kiss_size; i++) kiss_seed[i] = kiss_default_seed[i]; + } if (size != NULL) - *size = kiss_size; + *size = lehmer_lcg_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 (GFC_DESCRIPTOR_EXTENT(put,0) < kiss_size) - runtime_error ("Array size of PUT is too small."); - - /* We copy the seed given by the user. */ - for (i = 0; i < kiss_size; i++) - memcpy (seed + i * sizeof(GFC_UINTEGER_4), - &(put->base_addr[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(put,0)]), - sizeof(GFC_UINTEGER_4)); - - /* We put it after scrambling the bytes, to paper around users who - provide seeds with quality only in the lower or upper part. */ - scramble_seed ((unsigned char *) kiss_seed, seed, 4*kiss_size); + /* We copy the seed given by the user. */ + lehmer_lcg_seed = put->base_addr[0]; + if (lehmer_lcg_seed == 0) { + for (i = 0; i < kiss_size; i++) + kiss_seed[i] = kiss_default_seed[i]; + } else { + for (i = 0, seed = lehmer_lcg_seed; i < kiss_size; i++) { + seed = lehmer_lcg(seed); + kiss_seed[i] = seed; + } + } } /* 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 (GFC_DESCRIPTOR_EXTENT(get,0) < kiss_size) - runtime_error ("Array size of GET is too small."); - - /* Unscramble the seed. */ - unscramble_seed (seed, (unsigned char *) kiss_seed, 4*kiss_size); - - /* Then copy it back to the user variable. */ - for (i = 0; i < kiss_size; i++) - memcpy (&(get->base_addr[(kiss_size - 1 - i) * GFC_DESCRIPTOR_STRIDE(get,0)]), - seed + i * sizeof(GFC_UINTEGER_4), - sizeof(GFC_UINTEGER_4)); - } + get->base_addr[0] = lehmer_lcg_seed; __gthread_mutex_unlock (&random_lock); } @@ -734,6 +702,7 @@ void random_seed_i8 (GFC_INTEGER_8 *size, gfc_array_i8 *put, gfc_array_i8 *get) { int i; + GFC_UINTEGER_4 seed; __gthread_mutex_lock (&random_lock); @@ -743,45 +712,33 @@ random_seed_i8 (GFC_INTEGER_8 *size, gfc /* 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) + if (size == NULL && put == NULL && get == NULL) { + lehmer_lcg_seed = 0; for (i = 0; i < kiss_size; i++) kiss_seed[i] = kiss_default_seed[i]; + } if (size != NULL) - *size = kiss_size / 2; + *size = lehmer_lcg_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 (GFC_DESCRIPTOR_EXTENT(put,0) < 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 / 2; i++) - memcpy (&kiss_seed[2*i], &(put->base_addr[i * GFC_DESCRIPTOR_STRIDE(put,0)]), - sizeof (GFC_UINTEGER_8)); + /* We copy the seed given by the user. */ + lehmer_lcg_seed = put->base_addr[0]; + if (lehmer_lcg_seed == 0) { + for (i = 0; i < kiss_size; i++) + kiss_seed[i] = kiss_default_seed[i]; + } else { + for (i = 0, seed = lehmer_lcg_seed; i < kiss_size; i++) { + seed = lehmer_lcg(seed); + kiss_seed[i] = seed; + } + } } /* 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 (GFC_DESCRIPTOR_EXTENT(get,0) < 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 / 2; i++) - memcpy (&(get->base_addr[i * GFC_DESCRIPTOR_STRIDE(get,0)]), &kiss_seed[2*i], - sizeof (GFC_UINTEGER_8)); - } + get->base_addr[0] = lehmer_lcg_seed; __gthread_mutex_unlock (&random_lock); }
