src/z-rand.c

   1 /* File: z-rand.c */
   2
   3 /*
   4  * Copyright (c) 1997 Ben Harrison, and others
   5  *
   6  * This software may be copied and distributed for educational, research,
   7  * and not for profit purposes provided that this copyright and statement
   8  * are included in all such copies.  Other copyrights may also apply.
   9  */
  10
  11
  12 /*
  13  * This file provides an optimized random number generator.
  14  *
  15  *
  16  * This code provides both a "quick" random number generator (4 bytes of
  17  * state), and a "decent" random number generator (256 bytes of state),
  18  * both available in two flavors, first, the simple "mod" flavor, which
  19  * is fast, but slightly biased at high values, and second, the simple
  20  * "div" flavor, which is less fast (and potentially non-terminating)
  21  * but which is not biased and is much less subject to non-randomness
  22  * problems in the low bits.  Note the "rand_int()" macro in "z-rand.h",
  23  * which must specify a "default" flavor.
  24  *
  25  * Note the use of the "simple" RNG, first you activate it via
  26  * "Rand_quick = TRUE" and "Rand_value = seed" and then it is used
  27  * automatically used instead of the "complex" RNG, and when you are
  28  * done, you de-activate it via "Rand_quick = FALSE" or choose a new
  29  * seed via "Rand_value = seed".
  30  *
  31  *
  32  * This (optimized) random number generator is based loosely on the old
  33  * "random.c" file from Berkeley but with some major optimizations and
  34  * algorithm changes.  See below for more details.
  35  *
  36  * Some code by Ben Harrison (benh@phial.com).
  37  *
  38  * Some code by Randy (randy@stat.tamu.edu).
  39  */
  40
  41
  42
  43 #include "z-rand.h"
  44
  45
  46 /*
  47  * Random Number Generator -- Linear Congruent RNG
  48  */
  49 #define LCRNG(X)        ((X) * 1103515245 + 12345)
  50
  51
  52
  53 /*
  54  * Use the "simple" LCRNG
  55  */
  56 bool Rand_quick = TRUE;
  57
  58
  59 /*
  60  * Current "value" of the "simple" RNG
  61  */
  62 u32b Rand_value;
  63
  64
  65 /*
  66  * Current "index" for the "complex" RNG
  67  */
  68 u16b Rand_place;
  69
  70 /*
  71  * Current "state" table for the "complex" RNG
  72  */
  73 u32b Rand_state[RAND_DEG];
  74
  75
  76
  77 /*
  78  * Initialize the "complex" RNG using a new seed
  79  */
  80 void Rand_state_init(u32b seed)
  81 {
  82         int i, j;
  83
  84         /* Seed the table */
  85         Rand_state[0] = seed;
  86
  87         /* Propagate the seed */
  88         for (i = 1; i < RAND_DEG; i++) Rand_state[i] = LCRNG(Rand_state[i-1]);
  89
  90         /* Cycle the table ten times per degree */
  91         for (i = 0; i < RAND_DEG * 10; i++)
  92         {
  93                 /* Acquire the next index */
  94                 j = Rand_place + 1;
  95                 if (j == RAND_DEG) j = 0;
  96
  97                 /* Update the table, extract an entry */
  98                 Rand_state[j] += Rand_state[Rand_place];
  99
 100                 /* Advance the index */
 101                 Rand_place = j;
 102         }
 103 }
 104
 105
 106 /*
 107  * Extract a "random" number from 0 to m-1, via "division"
 108  *
 109  * This method selects "random" 28-bit numbers, and then uses
 110  * division to drop those numbers into "m" different partitions,
 111  * plus a small non-partition to reduce bias, taking as the final
 112  * value the first "good" partition that a number falls into.
 113  *
 114  * This method has no bias, and is much less affected by patterns
 115  * in the "low" bits of the underlying RNG's.
 116  *
 117  * Note that "m" must not be greater than 0x1000000, or division
 118  * by zero will result.
 119  *
 120  * ToDo: Check for m > 0x1000000.
 121  */
 122 u32b Rand_div(u32b m)
 123 {
 124         u32b r, n;
 125
 126         /* Hack -- simple case */
 127         if (m <= 1) return (0);
 128
 129         /* Partition size */
 130         n = (0x10000000 / m);
 131
 132         /* Use a simple RNG */
 133         if (Rand_quick)
 134         {
 135                 /* Wait for it */
 136                 while (1)
 137                 {
 138                         /* Cycle the generator */
 139                         r = (Rand_value = LCRNG(Rand_value));
 140
 141                         /* Mutate a 28-bit "random" number */
 142                         r = ((r >> 4) & 0x0FFFFFFF) / n;
 143
 144                         /* Done */
 145                         if (r < m) break;
 146                 }
 147         }
 148
 149         /* Use a complex RNG */
 150         else
 151         {
 152                 /* Wait for it */
 153                 while (1)
 154                 {
 155                         int j;
 156
 157                         /* Acquire the next index */
 158                         j = Rand_place + 1;
 159                         if (j == RAND_DEG) j = 0;
 160
 161                         /* Update the table, extract an entry */
 162                         r = (Rand_state[j] += Rand_state[Rand_place]);
 163
 164                         /* Hack -- extract a 28-bit "random" number */
 165                         r = ((r >> 4) & 0x0FFFFFFF) / n;
 166
 167                         /* Advance the index */
 168                         Rand_place = j;
 169
 170                         /* Done */
 171                         if (r < m) break;
 172                 }
 173         }
 174
 175         /* Use the value */
 176         return (r);
 177 }
 178
 179
 180
 181
 182 /*
 183  * The number of entries in the "Rand_normal_table"
 184  */
 185 #define RANDNOR_NUM     256
 186
 187 /*
 188  * The standard deviation of the "Rand_normal_table"
 189  */
 190 #define RANDNOR_STD     64
 191
 192 /*
 193  * The normal distribution table for the "Rand_normal()" function (below)
 194  */
 195 static s16b Rand_normal_table[RANDNOR_NUM] =
 196 {
 197         206,     613,    1022,    1430,         1838,    2245,    2652,    3058,
 198         3463,    3867,    4271,    4673,        5075,    5475,    5874,    6271,
 199         6667,    7061,    7454,    7845,        8234,    8621,    9006,    9389,
 200         9770,   10148,   10524,   10898,   11269,       11638,   12004,   12367,
 201         12727,   13085,   13440,   13792,   14140,      14486,   14828,   15168,
 202         15504,   15836,   16166,   16492,   16814,      17133,   17449,   17761,
 203         18069,   18374,   18675,   18972,   19266,      19556,   19842,   20124,
 204         20403,   20678,   20949,   21216,   21479,      21738,   21994,   22245,
 205
 206         22493,   22737,   22977,   23213,   23446,      23674,   23899,   24120,
 207         24336,   24550,   24759,   24965,   25166,      25365,   25559,   25750,
 208         25937,   26120,   26300,   26476,   26649,      26818,   26983,   27146,
 209         27304,   27460,   27612,   27760,   27906,      28048,   28187,   28323,
 210         28455,   28585,   28711,   28835,   28955,      29073,   29188,   29299,
 211         29409,   29515,   29619,   29720,   29818,      29914,   30007,   30098,
 212         30186,   30272,   30356,   30437,   30516,      30593,   30668,   30740,
 213         30810,   30879,   30945,   31010,   31072,      31133,   31192,   31249,
 214
 215         31304,   31358,   31410,   31460,   31509,      31556,   31601,   31646,
 216         31688,   31730,   31770,   31808,   31846,      31882,   31917,   31950,
 217         31983,   32014,   32044,   32074,   32102,      32129,   32155,   32180,
 218         32205,   32228,   32251,   32273,   32294,      32314,   32333,   32352,
 219         32370,   32387,   32404,   32420,   32435,      32450,   32464,   32477,
 220         32490,   32503,   32515,   32526,   32537,      32548,   32558,   32568,
 221         32577,   32586,   32595,   32603,   32611,      32618,   32625,   32632,
 222         32639,   32645,   32651,   32657,   32662,      32667,   32672,   32677,
 223
 224         32682,   32686,   32690,   32694,   32698,      32702,   32705,   32708,
 225         32711,   32714,   32717,   32720,   32722,      32725,   32727,   32729,
 226         32731,   32733,   32735,   32737,   32739,      32740,   32742,   32743,
 227         32745,   32746,   32747,   32748,   32749,      32750,   32751,   32752,
 228         32753,   32754,   32755,   32756,   32757,      32757,   32758,   32758,
 229         32759,   32760,   32760,   32761,   32761,      32761,   32762,   32762,
 230         32763,   32763,   32763,   32764,   32764,      32764,   32764,   32765,
 231         32765,   32765,   32765,   32766,   32766,      32766,   32766,   32767,
 232 };
 233
 234
 235
 236 /*
 237  * Generate a random integer number of NORMAL distribution
 238  *
 239  * The table above is used to generate a psuedo-normal distribution,
 240  * in a manner which is much faster than calling a transcendental
 241  * function to calculate a true normal distribution.
 242  *
 243  * Basically, entry 64*N in the table above represents the number of
 244  * times out of 32767 that a random variable with normal distribution
 245  * will fall within N standard deviations of the mean.  That is, about
 246  * 68 percent of the time for N=1 and 95 percent of the time for N=2.
 247  *
 248  * The table above contains a "faked" final entry which allows us to
 249  * pretend that all values in a normal distribution are strictly less
 250  * than four standard deviations away from the mean.  This results in
 251  * "conservative" distribution of approximately 1/32768 values.
 252  *
 253  * Note that the binary search takes up to 16 quick iterations.
 254  */
 255 s16b Rand_normal(int mean, int stand)
 256 {
 257         s16b tmp;
 258         s16b offset;
 259
 260         s16b low = 0;
 261         s16b high = RANDNOR_NUM;
 262
 263         /* Paranoia */
 264         if (stand < 1) return (mean);
 265
 266         /* Roll for probability */
 267         tmp = (s16b)rand_int(32768);
 268
 269         /* Binary Search */
 270         while (low < high)
 271         {
 272                 int mid = (low + high) >> 1;
 273
 274                 /* Move right if forced */
 275                 if (Rand_normal_table[mid] < tmp)
 276                 {
 277                         low = mid + 1;
 278                 }
 279
 280                 /* Move left otherwise */
 281                 else
 282                 {
 283                         high = mid;
 284                 }
 285         }
 286
 287         /* Convert the index into an offset */
 288         offset = (long)stand * (long)low / RANDNOR_STD;
 289
 290         /* One half should be negative */
 291         if (rand_int(100) < 50) return (mean - offset);
 292
 293         /* One half should be positive */
 294         return (mean + offset);
 295 }
 296
 297
 298 /*
 299  * Extract a "random" number from 0 to m-1, using the "simple" RNG.
 300  *
 301  * This function should be used when generating random numbers in
 302  * "external" program parts like the main-*.c files.  It preserves
 303  * the current RNG state to prevent influences on game-play.
 304  *
 305  * Could also use rand() from <stdlib.h> directly. XXX XXX XXX
 306  */
 307 u32b Rand_simple(u32b m)
 308 {
 309         static bool initialized = FALSE;
 310         static u32b simple_rand_value;
 311         bool old_rand_quick;
 312         u32b old_rand_value;
 313         u32b result;
 314
 315
 316         /* Save RNG state */
 317         old_rand_quick = Rand_quick;
 318         old_rand_value = Rand_value;
 319
 320         /* Use "simple" RNG */
 321         Rand_quick = TRUE;
 322
 323         if (initialized)
 324         {
 325                 /* Use stored seed */
 326                 Rand_value = simple_rand_value;
 327         }
 328         else
 329         {
 330                 /* Initialize with new seed */
 331                 Rand_value = time(NULL);
 332                 initialized = TRUE;
 333         }
 334
 335         /* Get a random number */
 336         result = rand_int(m);
 337
 338         /* Store the new seed */
 339         simple_rand_value = Rand_value;
 340
 341         /* Restore RNG state */
 342         Rand_quick = old_rand_quick;
 343         Rand_value = old_rand_value;
 344
 345         /* Use the value */
 346         return (result);
 347 }