utility/rand.c

   1 /**********************************************************************
   2  Freeciv - Copyright (C) 1996 - A Kjeldberg, L Gregersen, P Unold
   3    This program is free software; you can redistribute it and/or modify
   4    it under the terms of the GNU General Public License as published by
   5    the Free Software Foundation; either version 2, or (at your option)
   6    any later version.
   7
   8    This program is distributed in the hope that it will be useful,
   9    but WITHOUT ANY WARRANTY; without even the implied warranty of
  10    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11    GNU General Public License for more details.
  12 ***********************************************************************/
  13
  14 /*************************************************************************
  15    The following random number generator can be found in _The Art of
  16    Computer Programming Vol 2._ (2nd ed) by Donald E. Knuth. (C)  1998.
  17    The algorithm is described in section 3.2.2 as Mitchell and Moore's
  18    variant of a standard additive number generator.  Note that the
  19    the constants 55 and 24 are not random.  Please become familiar with
  20    this algorithm before you mess with it.
  21
  22    Since the additive number generator requires a table of numbers from
  23    which to generate its random sequences, we must invent a way to
  24    populate that table from a single seed value.  I have chosen to do
  25    this with a different PRNG, known as the "linear congruential method"
  26    (also found in Knuth, Vol2).  I must admit that my choices of constants
  27    (3, 257, and MAX_UINT32) are probably not optimal, but they seem to
  28    work well enough for our purposes.
  29
  30    Original author for this code: Cedric Tefft <cedric@earthling.net>
  31    Modified to use rand_state struct by David Pfitzner <dwp@mso.anu.edu.au>
  32 *************************************************************************/
  33
  34 #ifdef HAVE_CONFIG_H
  35 #include <fc_config.h>
  36 #endif
  37
  38 /* utility */
  39 #include "log.h"
  40 #include "support.h"            /* TRUE, FALSE */
  41
  42 #include "rand.h"
  43
  44 #define log_rand log_debug
  45
  46 /* A global random state:
  47  * Initialized by fc_srand(), updated by fc_rand(),
  48  * Can be duplicated/saved/restored via fc_rand_state()
  49  * and fc_rand_set_state().
  50  */
  51 static RANDOM_STATE rand_state;
  52
  53 /*************************************************************************
  54   Returns a new random value from the sequence, in the interval 0 to
  55   (size-1) inclusive, and updates global state for next call.
  56   This means that if size <= 1 the function will always return 0.
  57
  58   Once we calculate new_rand below uniform (we hope) between 0 and
  59   MAX_UINT32 inclusive, need to reduce to required range.  Using
  60   modulus is bad because generators like this are generally less
  61   random for their low-significance bits, so this can give poor
  62   results when 'size' is small.  Instead want to divide the range
  63   0..MAX_UINT32 into (size) blocks, each with (divisor) values, and
  64   for any remainder, repeat the calculation of new_rand.
  65   Then:
  66          return_val = new_rand / divisor;
  67   Will repeat for new_rand > max, where:
  68          max = size * divisor - 1
  69   Then max <= MAX_UINT32 implies
  70          size * divisor <= (MAX_UINT32+1)
  71   thus   divisor <= (MAX_UINT32+1)/size
  72
  73   Need to calculate this divisor.  Want divisor as large as possible
  74   given above contraint, but it doesn't hurt us too much if it is a
  75   bit smaller (just have to repeat more often).  Calculation exactly
  76   as above is complicated by fact that (MAX_UINT32+1) may not be
  77   directly representable in type RANDOM_TYPE, so we do instead:
  78          divisor = MAX_UINT32/size
  79 *************************************************************************/
  80 RANDOM_TYPE fc_rand_debug(RANDOM_TYPE size, const char *called_as,
  81                           int line, const char *file)
  82 {
  83   RANDOM_TYPE new_rand, divisor, max;
  84   int bailout = 0;
  85
  86   fc_assert_ret_val(rand_state.is_init, 0);
  87
  88   if (size > 1) {
  89     divisor = MAX_UINT32 / size;
  90     max = size * divisor - 1;
  91   } else {
  92     /* size == 0 || size == 1 */
  93
  94     /*
  95      * These assignments are only here to make the compiler
  96      * happy. Since each usage is guarded with a if(size>1).
  97      */
  98     max = MAX_UINT32;
  99     divisor = 1;
 100   }
 101
 102   do {
 103     new_rand = (rand_state.v[rand_state.j]
 104                 + rand_state.v[rand_state.k]) & MAX_UINT32;
 105
 106     rand_state.x = (rand_state.x +1) % 56;
 107     rand_state.j = (rand_state.j +1) % 56;
 108     rand_state.k = (rand_state.k +1) % 56;
 109     rand_state.v[rand_state.x] = new_rand;
 110
 111     if (++bailout > 10000) {
 112       log_error("%s(%lu) = %lu bailout at %s:%d",
 113                 called_as, (unsigned long) size,
 114                 (unsigned long) new_rand, file, line);
 115       new_rand = 0;
 116       break;
 117     }
 118
 119   } while (size > 1 && new_rand > max);
 120
 121   if (size > 1) {
 122     new_rand /= divisor;
 123   } else {
 124     new_rand = 0;
 125   }
 126
 127   log_rand("%s(%lu) = %lu at %s:%d",
 128            called_as, (unsigned long) size,
 129            (unsigned long) new_rand, file, line);
 130
 131   return new_rand;
 132 }
 133
 134 /*************************************************************************
 135   Initialize the generator; see comment at top of file.
 136 *************************************************************************/
 137 void fc_srand(RANDOM_TYPE seed)
 138 {
 139     int  i;
 140
 141     rand_state.v[0]=(seed & MAX_UINT32);
 142
 143     for(i=1; i<56; i++) {
 144        rand_state.v[i] = (3 * rand_state.v[i-1] + 257) & MAX_UINT32;
 145     }
 146
 147     rand_state.j = (55-55);
 148     rand_state.k = (55-24);
 149     rand_state.x = (55-0);
 150
 151     rand_state.is_init = TRUE;
 152
 153     /* Heat it up a bit:
 154      * Using modulus in fc_rand() this was important to pass
 155      * test_random1().  Now using divisor in fc_rand() that particular
 156      * test no longer indicates problems, but this seems a good idea
 157      * anyway -- eg, other tests could well reveal other initial
 158      * problems even using divisor.
 159      */
 160     for (i=0; i<10000; i++) {
 161       (void) fc_rand(MAX_UINT32);
 162     }
 163 }
 164
 165 /*************************************************************************
 166   Return whether the current state has been initialized.
 167 *************************************************************************/
 168 bool fc_rand_is_init(void)
 169 {
 170   return rand_state.is_init;
 171 }
 172
 173 /*************************************************************************
 174   Return a copy of the current rand_state; eg for save/restore.
 175 *************************************************************************/
 176 RANDOM_STATE fc_rand_state(void)
 177 {
 178   int i;
 179
 180   log_rand("fc_rand_state J=%d K=%d X=%d",
 181            rand_state.j, rand_state.k, rand_state.x);
 182   for (i  = 0; i < 8; i++) {
 183     log_rand("fc_rand_state %d, %08x %08x %08x %08x %08x %08x %08x",
 184              i, rand_state.v[7 * i],
 185              rand_state.v[7 * i + 1], rand_state.v[7 * i + 2],
 186              rand_state.v[7 * i + 3], rand_state.v[7 * i + 4],
 187              rand_state.v[7 * i + 5], rand_state.v[7 * i + 6]);
 188   }
 189
 190   return rand_state;
 191 }
 192
 193 /*************************************************************************
 194   Replace current rand_state with user-supplied; eg for save/restore.
 195   Caller should take care to set state.is_init beforehand if necessary.
 196 *************************************************************************/
 197 void fc_rand_set_state(RANDOM_STATE state)
 198 {
 199   int i;
 200
 201   rand_state = state;
 202
 203   log_rand("fc_rand_set_state J=%d K=%d X=%d",
 204            rand_state.j, rand_state.k, rand_state.x);
 205   for (i  = 0; i < 8; i++) {
 206     log_rand("fc_rand_set_state %d, %08x %08x %08x %08x %08x %08x %08x",
 207              i, rand_state.v[7 * i],
 208              rand_state.v[7 * i + 1], rand_state.v[7 * i + 2],
 209              rand_state.v[7 * i + 3], rand_state.v[7 * i + 4],
 210              rand_state.v[7 * i + 5], rand_state.v[7 * i + 6]);
 211   }
 212 }
 213
 214 /*************************************************************************
 215   Test one aspect of randomness, using n numbers.
 216   Reports results to LOG_TEST; with good randomness, behaviourchange
 217   and behavioursame should be about the same size.
 218   Tests current random state; saves and restores state, so can call
 219   without interrupting current sequence.
 220 *************************************************************************/
 221 void test_random1(int n)
 222 {
 223   RANDOM_STATE saved_state;
 224   int i, old_value = 0, new_value;
 225   bool didchange, olddidchange = FALSE;
 226   int behaviourchange = 0, behavioursame = 0;
 227
 228   saved_state = fc_rand_state();
 229   /* fc_srand(time(NULL)); */  /* use current state */
 230
 231   for (i = 0; i < n+2; i++) {
 232     new_value = fc_rand(2);
 233     if (i > 0) {                /* have old */
 234       didchange = (new_value != old_value);
 235       if (i > 1) {              /* have olddidchange */
 236         if (didchange != olddidchange) {
 237           behaviourchange++;
 238         } else {
 239           behavioursame++;
 240         }
 241       }
 242       olddidchange = didchange;
 243     }
 244     old_value = new_value;
 245   }
 246   log_test("test_random1(%d) same: %d, change: %d",
 247            n, behavioursame, behaviourchange);
 248
 249   /* restore state: */
 250   fc_rand_set_state(saved_state);
 251 }
 252
 253 /*************************************************************************
 254   Local pseudo-random function for repeatedly reaching the same result,
 255   instead of fc_rand().  Primarily needed for tiles.
 256
 257   Use an invariant equation for seed.
 258   Result is 0 to (size - 1).
 259 *************************************************************************/
 260 RANDOM_TYPE fc_randomly_debug(RANDOM_TYPE seed, RANDOM_TYPE size,
 261                               const char *called_as,
 262                               int line, const char *file)
 263 {
 264   RANDOM_TYPE result;
 265
 266 #define LARGE_PRIME (10007)
 267 #define SMALL_PRIME (1009)
 268
 269   /* Check for overflow and underflow */
 270   fc_assert_ret_val(seed < MAX_UINT32 / LARGE_PRIME, 0);
 271   fc_assert_ret_val(size < SMALL_PRIME, 0);
 272   fc_assert_ret_val(size > 0, 0);
 273   result = ((seed * LARGE_PRIME) % SMALL_PRIME) % size;
 274
 275   log_rand("%s(%lu,%lu) = %lu at %s:%d",
 276            called_as, (unsigned long) seed, (unsigned long) size,
 277            (unsigned long) result, file, line);
 278
 279   return result;
 280 }