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 "shared.h"
  41 #include "support.h"            /* TRUE, FALSE */
  42
  43 #include "rand.h"
  44
  45 #define log_rand log_debug
  46
  47 /* A global random state:
  48  * Initialized by fc_srand(), updated by fc_rand(),
  49  * Can be duplicated/saved/restored via fc_rand_state()
  50  * and fc_rand_set_state().
  51  */
  52 static RANDOM_STATE rand_state;
  53
  54 /*************************************************************************
  55   Returns a new random value from the sequence, in the interval 0 to
  56   (size-1) inclusive, and updates global state for next call.
  57   This means that if size <= 1 the function will always return 0.
  58
  59   Once we calculate new_rand below uniform (we hope) between 0 and
  60   MAX_UINT32 inclusive, need to reduce to required range.  Using
  61   modulus is bad because generators like this are generally less
  62   random for their low-significance bits, so this can give poor
  63   results when 'size' is small.  Instead want to divide the range
  64   0..MAX_UINT32 into (size) blocks, each with (divisor) values, and
  65   for any remainder, repeat the calculation of new_rand.
  66   Then:
  67          return_val = new_rand / divisor;
  68   Will repeat for new_rand > max, where:
  69          max = size * divisor - 1
  70   Then max <= MAX_UINT32 implies
  71          size * divisor <= (MAX_UINT32+1)
  72   thus   divisor <= (MAX_UINT32+1)/size
  73
  74   Need to calculate this divisor.  Want divisor as large as possible
  75   given above contraint, but it doesn't hurt us too much if it is a
  76   bit smaller (just have to repeat more often).  Calculation exactly
  77   as above is complicated by fact that (MAX_UINT32+1) may not be
  78   directly representable in type RANDOM_TYPE, so we do instead:
  79          divisor = MAX_UINT32/size
  80 *************************************************************************/
  81 RANDOM_TYPE fc_rand_debug(RANDOM_TYPE size, const char *called_as,
  82                           int line, const char *file)
  83 {
  84   RANDOM_TYPE new_rand, divisor, max;
  85   int bailout = 0;
  86
  87   fc_assert_ret_val(rand_state.is_init, 0);
  88
  89   if (size > 1) {
  90     divisor = MAX_UINT32 / size;
  91     max = size * divisor - 1;
  92   } else {
  93     /* size == 0 || size == 1 */
  94
  95     /*
  96      * These assignments are only here to make the compiler
  97      * happy. Since each usage is guarded with a if (size > 1).
  98      */
  99     max = MAX_UINT32;
 100     divisor = 1;
 101   }
 102
 103   do {
 104     new_rand = (rand_state.v[rand_state.j]
 105                 + rand_state.v[rand_state.k]) & MAX_UINT32;
 106
 107     rand_state.x = (rand_state.x +1) % 56;
 108     rand_state.j = (rand_state.j +1) % 56;
 109     rand_state.k = (rand_state.k +1) % 56;
 110     rand_state.v[rand_state.x] = new_rand;
 111
 112     if (++bailout > 10000) {
 113       log_error("%s(%lu) = %lu bailout at %s:%d",
 114                 called_as, (unsigned long) size,
 115                 (unsigned long) new_rand, file, line);
 116       new_rand = 0;
 117       break;
 118     }
 119
 120   } while (size > 1 && new_rand > max);
 121
 122   if (size > 1) {
 123     new_rand /= divisor;
 124   } else {
 125     new_rand = 0;
 126   }
 127
 128   log_rand("%s(%lu) = %lu at %s:%d",
 129            called_as, (unsigned long) size,
 130            (unsigned long) new_rand, file, line);
 131
 132   return new_rand;
 133 }
 134
 135 /*************************************************************************
 136   Initialize the generator; see comment at top of file.
 137 *************************************************************************/
 138 void fc_srand(RANDOM_TYPE seed)
 139 {
 140   int  i;
 141
 142   rand_state.v[0] = (seed & MAX_UINT32);
 143
 144   for (i = 1; i < 56; i++) {
 145     rand_state.v[i] = (3 * rand_state.v[i-1] + 257) & MAX_UINT32;
 146   }
 147
 148   rand_state.j = (55-55);
 149   rand_state.k = (55-24);
 150   rand_state.x = (55-0);
 151
 152   rand_state.is_init = TRUE;
 153
 154   /* Heat it up a bit:
 155    * Using modulus in fc_rand() this was important to pass
 156    * test_random1().  Now using divisor in fc_rand() that particular
 157    * test no longer indicates problems, but this seems a good idea
 158    * anyway -- eg, other tests could well reveal other initial
 159    * problems even using divisor.
 160    */
 161   for (i = 0; i < 10000; i++) {
 162     (void) fc_rand(MAX_UINT32);
 163   }
 164 }
 165
 166 /*************************************************************************
 167   Return whether the current state has been initialized.
 168 *************************************************************************/
 169 bool fc_rand_is_init(void)
 170 {
 171   return rand_state.is_init;
 172 }
 173
 174 /*************************************************************************
 175   Return a copy of the current rand_state; eg for save/restore.
 176 *************************************************************************/
 177 RANDOM_STATE fc_rand_state(void)
 178 {
 179   int i;
 180
 181   log_rand("fc_rand_state J=%d K=%d X=%d",
 182            rand_state.j, rand_state.k, rand_state.x);
 183   for (i  = 0; i < 8; i++) {
 184     log_rand("fc_rand_state %d, %08x %08x %08x %08x %08x %08x %08x",
 185              i, rand_state.v[7 * i],
 186              rand_state.v[7 * i + 1], rand_state.v[7 * i + 2],
 187              rand_state.v[7 * i + 3], rand_state.v[7 * i + 4],
 188              rand_state.v[7 * i + 5], rand_state.v[7 * i + 6]);
 189   }
 190
 191   return rand_state;
 192 }
 193
 194 /*************************************************************************
 195   Replace current rand_state with user-supplied; eg for save/restore.
 196   Caller should take care to set state.is_init beforehand if necessary.
 197 *************************************************************************/
 198 void fc_rand_set_state(RANDOM_STATE state)
 199 {
 200   int i;
 201
 202   rand_state = state;
 203
 204   log_rand("fc_rand_set_state J=%d K=%d X=%d",
 205            rand_state.j, rand_state.k, rand_state.x);
 206   for (i  = 0; i < 8; i++) {
 207     log_rand("fc_rand_set_state %d, %08x %08x %08x %08x %08x %08x %08x",
 208              i, rand_state.v[7 * i],
 209              rand_state.v[7 * i + 1], rand_state.v[7 * i + 2],
 210              rand_state.v[7 * i + 3], rand_state.v[7 * i + 4],
 211              rand_state.v[7 * i + 5], rand_state.v[7 * i + 6]);
 212   }
 213 }
 214
 215 /*************************************************************************
 216   Test one aspect of randomness, using n numbers.
 217   Reports results to LOG_TEST; with good randomness, behaviourchange
 218   and behavioursame should be about the same size.
 219   Tests current random state; saves and restores state, so can call
 220   without interrupting current sequence.
 221 *************************************************************************/
 222 void test_random1(int n)
 223 {
 224   RANDOM_STATE saved_state;
 225   int i, old_value = 0, new_value;
 226   bool didchange, olddidchange = FALSE;
 227   int behaviourchange = 0, behavioursame = 0;
 228
 229   saved_state = fc_rand_state();
 230   /* fc_srand(time(NULL)); */  /* use current state */
 231
 232   for (i = 0; i < n+2; i++) {
 233     new_value = fc_rand(2);
 234     if (i > 0) {                /* have old */
 235       didchange = (new_value != old_value);
 236       if (i > 1) {              /* have olddidchange */
 237         if (didchange != olddidchange) {
 238           behaviourchange++;
 239         } else {
 240           behavioursame++;
 241         }
 242       }
 243       olddidchange = didchange;
 244     }
 245     old_value = new_value;
 246   }
 247   log_test("test_random1(%d) same: %d, change: %d",
 248            n, behavioursame, behaviourchange);
 249
 250   /* restore state: */
 251   fc_rand_set_state(saved_state);
 252 }
 253
 254 /*************************************************************************
 255   Local pseudo-random function for repeatedly reaching the same result,
 256   instead of fc_rand().  Primarily needed for tiles.
 257
 258   Use an invariant equation for seed.
 259   Result is 0 to (size - 1).
 260 *************************************************************************/
 261 RANDOM_TYPE fc_randomly_debug(RANDOM_TYPE seed, RANDOM_TYPE size,
 262                               const char *called_as,
 263                               int line, const char *file)
 264 {
 265   RANDOM_TYPE result;
 266
 267 #define LARGE_PRIME (10007)
 268 #define SMALL_PRIME (1009)
 269
 270   /* Check for overflow and underflow */
 271   fc_assert_ret_val(seed < MAX_UINT32 / LARGE_PRIME, 0);
 272   fc_assert_ret_val(size < SMALL_PRIME, 0);
 273   fc_assert_ret_val(size > 0, 0);
 274   result = ((seed * LARGE_PRIME) % SMALL_PRIME) % size;
 275
 276   log_rand("%s(%lu,%lu) = %lu at %s:%d",
 277            called_as, (unsigned long) seed, (unsigned long) size,
 278            (unsigned long) result, file, line);
 279
 280   return result;
 281 }