[libgfortran,patch] Scramble bytes in the user-provided random seeds (PR 32812)

[FX Coudert <fxcoudert@gmail.com>]

[-- Attachment #1: Type: text/plain, Size: 1263 bytes --]

Hi all,

Current library implementation of RANDOM_NUMBER() is such that the  
highest bits of the user-provided seed determine the mots significant  
bits of floating-points values return, and the lowest bits of the  
seed control the least significant bits of fp values. Thus, a user  
providing us with a seed with good lower bits but poor higher bits  
will find that the fp values are not very random at all... one such  
example being the use of DATE_AND_TIME to provide for random seeds.

This patch adds a little to the quality of implementation by adding a  
(bijective) scrambling operation between the seed the user gives us  
and the seed we actually use. Two functions, scramble_seed and  
unscramble_seed, are used when  RANDOM_NUMBER() is called with a PUT=  
or GET= argument, which shuffles bytes of the seed.

Regtested on x86_64-linux, with both -m32 and -m64. Tested manually  
with different seeds of different qualities.


FX

PS: Of course, we're not adding quality to poor seeds (I'm not  
changing entropy with a bijective function!). We're just providing a  
little invisible help to one particular case of half-bad, half-good  
seeds.


-- 
François-Xavier Coudert
http://www.homepages.ucl.ac.uk/~uccafco/

[-- Attachment #2: random.ChangeLog --]
[-- Type: application/octet-stream, Size: 280 bytes --]

2008-02-25  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>

	PR libfortran/32812
	* intrinsics/random.c (scramble_seed, unscramble_seed): New
	functions.
	(random_seed_i4): Scramble the seed the user gives us before
	storing it, and unscramble it when we return it back later.


[-- Attachment #3: random.diff --]
[-- Type: application/octet-stream, Size: 2540 bytes --]

Index: intrinsics/random.c
===================================================================
--- intrinsics/random.c	(revision 132592)
+++ intrinsics/random.c	(working copy)
@@ -639,6 +639,29 @@ 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.  */
@@ -647,6 +670,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];
 
   __gthread_mutex_lock (&random_lock);
 
@@ -673,9 +697,15 @@ random_seed_i4 (GFC_INTEGER_4 *size, gfc
       if (((put->dim[0].ubound + 1 - put->dim[0].lbound)) < kiss_size)
         runtime_error ("Array size of PUT is too small.");
 
-      /*  This code now should do correct strides.  */
+      /*  We copy the seed given by the user.  */
       for (i = 0; i < kiss_size; i++)
-	kiss_seed[i] = (GFC_UINTEGER_4) put->data[i * put->dim[0].stride];
+	memcpy (seed + i * sizeof(GFC_UINTEGER_4),
+		&(put->data[(kiss_size - 1 - i) * put->dim[0].stride]),
+		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);
     }
 
   /* Return the seed to GET data.  */
@@ -689,9 +719,14 @@ random_seed_i4 (GFC_INTEGER_4 *size, gfc
       if (((get->dim[0].ubound + 1 - get->dim[0].lbound)) < kiss_size)
 	runtime_error ("Array size of GET is too small.");
 
-      /*  This code now should do correct strides.  */
+      /* 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++)
-        get->data[i * get->dim[0].stride] = (GFC_INTEGER_4) kiss_seed[i];
+       memcpy (&(get->data[(kiss_size - 1 - i) * get->dim[0].stride]),
+               seed + i * sizeof(GFC_UINTEGER_4),
+               sizeof(GFC_UINTEGER_4));
     }
 
   __gthread_mutex_unlock (&random_lock);