[PATCH, fortran/52879] RANDOM_SEED revisited
Steve Kargl
sgk@troutmask.apl.washington.edu
Sun Feb 9 00:40:00 GMT 2014
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);
}
More information about the Fortran
mailing list