[Freeciv-Dev] (PR#4417) Use of ISAAC Random generator in Freeciv

["Rafa³ Bursig" <bursig@xxxxxxxxx>]

Hi all

Here is implementation of ISAAC random number generator for freeciv.
It is faster from current code ( ~3x ) and give better (IMHO) result.
This patch is only proposition and don't replace old code fully 
(savegame code is needed)

Rafal

diff -u -r fc/common/rand.c fc_new/common/rand.c
--- fc/common/rand.c    Wed Jan 29 13:10:17 2003
+++ fc_new/common/rand.c        Sun Jun 15 22:09:43 2003
@@ -11,25 +11,14 @@
    GNU General Public License for more details.
 ***********************************************************************/
 
-/*************************************************************************
-   The following random number generator can be found in _The Art of 
-   Computer Programming Vol 2._ (2nd ed) by Donald E. Knuth. (C)  1998.
-   The algorithm is described in section 3.2.2 as Mitchell and Moore's
-   variant of a standard additive number generator.  Note that the
-   the constants 55 and 24 are not random.  Please become familiar with
-   this algorithm before you mess with it.
-
-   Since the additive number generator requires a table of numbers from
-   which to generate its random sequences, we must invent a way to 
-   populate that table from a single seed value.  I have chosen to do
-   this with a different PRNG, known as the "linear congruential method" 
-   (also found in Knuth, Vol2).  I must admit that my choices of constants
-   (3, 257, and MAX_UINT32) are probably not optimal, but they seem to
-   work well enough for our purposes.
-   
-   Original author for this code: Cedric Tefft <cedric@xxxxxxxxxxxxx>
-   Modified to use rand_state struct by David Pfitzner <dwp@xxxxxxxxxxxxxx>
-*************************************************************************/
+/*
+------------------------------------------------------------------------------
+ISAAC random number generator by Bob Jenkins.    Public Domain.
+http://burtleburtle.net/bob/rand/isaacafa.html
+
+this code may not work on 64 bit platforms -- Rafal
+------------------------------------------------------------------------------
+*/
 
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -50,111 +39,129 @@
  */
 static RANDOM_STATE rand_state;
 
-/*************************************************************************
-  Returns a new random value from the sequence, in the interval 0 to
-  (size-1) inclusive, and updates global state for next call.
-  This means that if size <= 1 the function will always return 0.
-
-  Once we calculate new_rand below uniform (we hope) between 0 and
-  MAX_UINT32 inclusive, need to reduce to required range.  Using
-  modulus is bad because generators like this are generally less
-  random for their low-significance bits, so this can give poor
-  results when 'size' is small.  Instead want to divide the range
-  0..MAX_UINT32 into (size) blocks, each with (divisor) values, and
-  for any remainder, repeat the calculation of new_rand.
-  Then:
-        return_val = new_rand / divisor;
-  Will repeat for new_rand > max, where:
-         max = size * divisor - 1
-  Then max <= MAX_UINT32 implies
-        size * divisor <= (MAX_UINT32+1)
-  thus   divisor <= (MAX_UINT32+1)/size
+
+#define ind(mm,x)  (*(RANDOM_TYPE *)((unsigned char *)(mm) + ((x) & 
((RANDSIZ-1)<<2))))
+
+#define rngstep(mix,a,b,mm,m,m2,r,x)   \
+do {                                   \
+  x = *m;                              \
+  a = (a^(mix)) + *(m2++);             \
+  *(m++) = y = ind(mm,x) + a + b;      \
+  *(r++) = b = ind(mm,y>>RANDSIZL) + x; \
+} while(0)
+
+static inline void isaac_rand(void)
+{
+   register RANDOM_TYPE a,b,x,y,*m,*mm,*m2,*r,*mend;
   
-  Need to calculate this divisor.  Want divisor as large as possible
-  given above contraint, but it doesn't hurt us too much if it is a
-  bit smaller (just have to repeat more often).  Calculation exactly
-  as above is complicated by fact that (MAX_UINT32+1) may not be
-  directly representable in type RANDOM_TYPE, so we do instead:
-         divisor = MAX_UINT32/size
+   mm = rand_state.randmem;
+   r  = rand_state.randrsl;
+   a  = rand_state.randa;
+   b  = rand_state.randb + (++(rand_state.randc));
+  
+   for (m = mm, mend = m2 = m+(RANDSIZ/2); m<mend; )
+   {
+      rngstep(a<<13, a, b, mm, m, m2, r, x);
+      rngstep(a>>6 , a, b, mm, m, m2, r, x);
+      rngstep(a<<2 , a, b, mm, m, m2, r, x);
+      rngstep(a>>16, a, b, mm, m, m2, r, x);
+   }
+   for (m2 = mm; m2<mend; )
+   {
+      rngstep(a<<13, a, b, mm, m, m2, r, x);
+      rngstep(a>>6 , a, b, mm, m, m2, r, x);
+      rngstep(a<<2 , a, b, mm, m, m2, r, x);
+      rngstep(a>>16, a, b, mm, m, m2, r, x);
+   }
+   
+   rand_state.randb = b;
+   rand_state.randa = a;
+}
+
+
+/*************************************************************************
+  Returns a new random value from the sequence.
 *************************************************************************/
 RANDOM_TYPE myrand(RANDOM_TYPE size) 
 {
-  RANDOM_TYPE new_rand, divisor, max;
-  int bailout = 0;
-
   assert(rand_state.is_init);
-    
-  if (size > 1) {
-    divisor = MAX_UINT32 / size;
-    max = size * divisor - 1;
-  } else {
-    /* size == 0 || size == 1 */
-
-    /* 
-     * These assignments are only here to make the compiler
-     * happy. Since each usage is guarded with a if(size>1).
-     */
-    max = MAX_UINT32;
-    divisor = 1;
-  }
-
-  do {
-    new_rand = (rand_state.v[rand_state.j]
-               + rand_state.v[rand_state.k]) & MAX_UINT32;
-
-    rand_state.x = (rand_state.x +1) % 56;
-    rand_state.j = (rand_state.j +1) % 56;
-    rand_state.k = (rand_state.k +1) % 56;
-    rand_state.v[rand_state.x] = new_rand;
-
-    if (++bailout > 10000) {
-      freelog(LOG_ERROR, "Bailout in myrand(%u)", size);
-      new_rand = 0;
-      break;
-    }
-
-  } while (size > 1 && new_rand > max);
-
-  if (size > 1) {
-    new_rand /= divisor;
-  } else {
-    new_rand = 0;
-  }
+  return ((!rand_state.randcnt-- ?
+     (isaac_rand(), rand_state.randcnt=RANDSIZ-1, 
rand_state.randrsl[rand_state.randcnt]) : \
+     rand_state.randrsl[rand_state.randcnt]) % (size + 1));
+} 
 
-  /* freelog(LOG_DEBUG, "rand(%u) = %u", size, new_rand); */
+#define mix(a,b,c,d,e,f,g,h)   \
+do {                           \
+   a^=b<<11; d+=a; b+=c;       \
+   b^=c>>2;  e+=b; c+=d;       \
+   c^=d<<8;  f+=c; d+=e;       \
+   d^=e>>16; g+=d; e+=f;       \
+   e^=f<<10; h+=e; f+=g;       \
+   f^=g>>4;  a+=f; g+=h;       \
+   g^=h<<8;  b+=g; h+=a;       \
+   h^=a>>9;  c+=h; a+=b;       \
+} while(0)
 
-  return new_rand;
-} 
 
 /*************************************************************************
   Initialize the generator; see comment at top of file.
 *************************************************************************/
 void mysrand(RANDOM_TYPE seed) 
 { 
-    int  i; 
-
-    rand_state.v[0]=(seed & MAX_UINT32);
-
-    for(i=1; i<56; i++) {
-       rand_state.v[i] = (3 * rand_state.v[i-1] + 257) & MAX_UINT32;
+    int i; 
+    RANDOM_TYPE a,b,c,d,e,f,g,h;
+    RANDOM_TYPE *m,*r;
+  
+    rand_state.randa = rand_state.randb = rand_state.randc = 0;
+    m = rand_state.randmem;
+    r = rand_state.randrsl;
+    a=b=c=d=e=f=g=h=0x9e3779b9;  /* the golden ratio */
+    
+    for (i=0; i<4; ++i)          /* scramble it */
+    {
+      mix(a,b,c,d,e,f,g,h);
     }
 
-    rand_state.j = (55-55);
-    rand_state.k = (55-24);
-    rand_state.x = (55-0);
+    if (seed) 
+    {
+      for (i=0; i<RANDSIZ; ++i) r[i] = seed;
+      for (i=0; i<RANDSIZ; i+=8)
+      {
+       
+
+         
+       /* initialize using the contents of r[] as the seed */  
+        a+=r[i  ]; b+=r[i+1]; c+=r[i+2]; d+=r[i+3];
+        e+=r[i+4]; f+=r[i+5]; g+=r[i+6]; h+=r[i+7];
+        mix(a,b,c,d,e,f,g,h);
+        m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
+        m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
+      }
+      /* do a second pass to make all of the seed affect all of m */
+      for (i=0; i<RANDSIZ; i+=8)
+      {
+        a+=m[i  ]; b+=m[i+1]; c+=m[i+2]; d+=m[i+3];
+        e+=m[i+4]; f+=m[i+5]; g+=m[i+6]; h+=m[i+7];
+        mix(a,b,c,d,e,f,g,h);
+        m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
+        m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
+      }
+    }
+    else
+    {
+      /* fill in mm[] with messy stuff */
+      for (i=0; i<RANDSIZ; i+=8)
+      {
+        mix(a,b,c,d,e,f,g,h);
+        m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
+        m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
+      }
+    }
 
+    isaac_rand();               /* fill in the first set of results */
+    rand_state.randcnt=RANDSIZ; /* prepare to use the first set of results */
     rand_state.is_init = TRUE;
-
-    /* Heat it up a bit:
-     * Using modulus in myrand() this was important to pass
-     * test_random1().  Now using divisor in myrand() that particular
-     * test no longer indicates problems, but this seems a good idea
-     * anyway -- eg, other tests could well reveal other initial
-     * problems even using divisor.
-     */
-    for (i=0; i<10000; i++) {
-      (void) myrand(MAX_UINT32);
-    }
+ 
 } 
 
 /*************************************************************************
@@ -197,8 +204,7 @@
   int behaviourchange = 0, behavioursame = 0;
 
   saved_state = get_myrand_state();
-  /* mysrand(time(NULL)); */  /* use current state */
-
+  
   for (i = 0; i < n+2; i++) {
     new_value = myrand(2);
     if (i > 0) {               /* have old */
diff -u -r fc/common/rand.h fc_new/common/rand.h
--- fc/common/rand.h    Fri Apr 12 15:50:56 2002
+++ fc_new/common/rand.h        Sun Jun 15 21:09:44 2003
@@ -20,9 +20,16 @@
 
 typedef unsigned int RANDOM_TYPE;
 
+#define RANDSIZL   (4)  /* I recommend 8 for crypto, 4 for simulations */
+#define RANDSIZ    (1<<RANDSIZL)
+
 typedef struct {
-  RANDOM_TYPE v[56];
-  int j, k, x;
+  RANDOM_TYPE randcnt;
+  RANDOM_TYPE randrsl[RANDSIZ];
+  RANDOM_TYPE randmem[RANDSIZ];
+  RANDOM_TYPE randa;
+  RANDOM_TYPE randb;
+  RANDOM_TYPE randc;
   bool is_init;                        /* initially 0 for static storage */
 } RANDOM_STATE;