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 /* Purpose: a simple random number generator -BEN- */
  13
  14 #if defined(WINDOWS)
  15 #include <Windows.h>
  16 #endif
  17
  18 #include "z-rand.h"
  19
  20
  21
  22
  23 /*
  24  * Angband 2.7.9 introduced a new (optimized) random number generator,
  25  * based loosely on the old "random.c" from Berkeley but with some major
  26  * optimizations and algorithm changes.  See below for more details.
  27  *
  28  * Code by myself (benh@phial.com) and Randy (randy@stat.tamu.edu).
  29  *
  30  * This code provides (1) a "decent" RNG, based on the "BSD-degree-63-RNG"
  31  * used in Angband 2.7.8, but rather optimized, and (2) a "simple" RNG,
  32  * based on the simple "LCRNG" currently used in Angband, but "corrected"
  33  * to give slightly better values.  Both of these are available in two
  34  * flavors, first, the simple "mod" flavor, which is fast, but slightly
  35  * biased at high values, and second, the simple "div" flavor, which is
  36  * less fast (and potentially non-terminating) but which is not biased
  37  * and is much less subject to low-bit-non-randomness problems.
  38  *
  39  * You can select your favorite flavor by proper definition of the
  40  * "randint0()" macro in the "defines.h" file.
  41  *
  42  * Note that, in Angband 2.8.0, the "state" table will be saved in the
  43  * savefile, so a special "initialization" phase will be necessary.
  44  *
  45  * Note the use of the "simple" RNG, first you activate it via
  46  * "Rand_quick = TRUE" and "Rand_value = seed" and then it is used
  47  * automatically used instead of the "complex" RNG, and when you are
  48  * done, you de-activate it via "Rand_quick = FALSE" or choose a new
  49  * seed via "Rand_value = seed".
  50  *
  51  *
  52  * RNG algorithm was fully rewritten. Upper comment is OLD.
  53  */
  54
  55
  56 /*
  57  * Currently unused
  58  */
  59 u16b Rand_place;
  60
  61 /*
  62  * Current "state" table for the RNG
  63  * Only index 0 to 3 are used
  64  */
  65 u32b Rand_state[RAND_DEG] = {
  66         123456789,
  67         362436069,
  68         521288629,
  69         88675123,
  70 };
  71
  72
  73 static u32b u32b_rotl(const u32b x, int k)
  74 {
  75         return (x << k) | (x >> (32 - k));
  76 }
  77
  78 /*
  79  * Initialize RNG state
  80  */
  81 static void Rand_seed(u32b seed, u32b* state)
  82 {
  83         int i;
  84
  85         for (i = 1; i <= 4; ++ i) {
  86                 seed = 1812433253UL * (seed ^ (seed >> 30)) + i;
  87                 state[i-1] = seed;
  88         }
  89 }
  90
  91 /*
  92  * Xoshiro128** Algorithm
  93  */
  94 static u32b Rand_Xoshiro128starstar(u32b *state)
  95 {
  96         const u32b result = u32b_rotl(state[0] * 5, 7) * 9;
  97
  98         const u32b t = state[1] << 9;
  99
 100         state[2] ^= state[0];
 101         state[3] ^= state[1];
 102         state[1] ^= state[2];
 103         state[0] ^= state[3];
 104
 105         state[2] ^= t;
 106
 107         state[3] = u32b_rotl(state[3], 11);
 108
 109         return result;
 110 }
 111
 112 static const u32b Rand_Xorshift_max = 0xFFFFFFFF;
 113
 114 /*
 115  * Initialize the RNG using a new seed
 116  */
 117 void Rand_state_set(u32b seed)
 118 {
 119         Rand_seed(seed, Rand_state);
 120 }
 121
 122 void Rand_state_init(void)
 123 {
 124 #ifdef RNG_DEVICE
 125
 126         FILE *fp = fopen(RNG_DEVICE, "r");
 127         int n;
 128
 129         do {
 130                 n = fread(Rand_state, sizeof(Rand_state[0]), 4, fp);
 131         } while (n != 4 || (Rand_state[0] | Rand_state[1] | Rand_state[2] | Rand_state[3]) == 0);
 132
 133         fclose(fp);
 134
 135 #elif defined(WINDOWS)
 136
 137         HCRYPTPROV hProvider;
 138
 139         CryptAcquireContext(&hProvider, NULL, NULL, PROV_RSA_FULL, 0);
 140
 141         do {
 142                 CryptGenRandom(hProvider, sizeof(Rand_state[0]) * 4, (BYTE*)Rand_state);
 143         } while ((Rand_state[0] | Rand_state[1] | Rand_state[2] | Rand_state[3]) == 0);
 144
 145         CryptReleaseContext(hProvider, 0);
 146
 147 #else
 148
 149         /* Basic seed */
 150         u32b seed = (time(NULL));
 151 #ifdef SET_UID
 152         /* Mutate the seed on Unix machines */
 153         seed = ((seed >> 3) * (getpid() << 1));
 154 #endif
 155         /* Seed the RNG */
 156         Rand_state_set(seed);
 157
 158 #endif
 159 }
 160
 161 /*
 162  * Backup the RNG state
 163  */
 164 void Rand_state_backup(u32b* backup_state)
 165 {
 166         int i;
 167
 168         for (i = 0; i < 4; ++ i) {
 169                 backup_state[i] = Rand_state[i];
 170         }
 171 }
 172
 173 /*
 174  * Restore the RNG state
 175  */
 176 void Rand_state_restore(u32b* backup_state)
 177 {
 178         int i;
 179
 180         for (i = 0; i < 4; ++ i) {
 181                 Rand_state[i] = backup_state[i];
 182         }
 183 }
 184
 185
 186 /*
 187  * Extract a "random" number from 0 to m-1, via "division"
 188  */
 189 static s32b Rand_div_impl(s32b m, u32b* state)
 190 {
 191         u32b scaling;
 192         u32b past;
 193         u32b ret;
 194
 195         /* Hack -- simple case */
 196         if (m <= 1) return (0);
 197
 198         scaling = Rand_Xorshift_max / m;
 199         past = scaling * m;
 200
 201         do {
 202                 ret = Rand_Xoshiro128starstar(state);
 203         } while (ret >= past);
 204
 205         return ret / scaling;
 206 }
 207
 208 s32b Rand_div(s32b m)
 209 {
 210         return Rand_div_impl(m, Rand_state);
 211 }
 212
 213
 214
 215
 216 /*
 217  * The number of entries in the "randnor_table"
 218  */
 219 #define RANDNOR_NUM     256
 220
 221 /*
 222  * The standard deviation of the "randnor_table"
 223  */
 224 #define RANDNOR_STD     64
 225
 226 /*
 227  * The normal distribution table for the "randnor()" function (below)
 228  */
 229 static s16b randnor_table[RANDNOR_NUM] =
 230 {
 231         206,     613,    1022,    1430,         1838,    2245,    2652,    3058,
 232         3463,    3867,    4271,    4673,        5075,    5475,    5874,    6271,
 233         6667,    7061,    7454,    7845,        8234,    8621,    9006,    9389,
 234         9770,   10148,   10524,   10898,   11269,       11638,   12004,   12367,
 235         12727,   13085,   13440,   13792,   14140,      14486,   14828,   15168,
 236         15504,   15836,   16166,   16492,   16814,      17133,   17449,   17761,
 237         18069,   18374,   18675,   18972,   19266,      19556,   19842,   20124,
 238         20403,   20678,   20949,   21216,   21479,      21738,   21994,   22245,
 239
 240         22493,   22737,   22977,   23213,   23446,      23674,   23899,   24120,
 241         24336,   24550,   24759,   24965,   25166,      25365,   25559,   25750,
 242         25937,   26120,   26300,   26476,   26649,      26818,   26983,   27146,
 243         27304,   27460,   27612,   27760,   27906,      28048,   28187,   28323,
 244         28455,   28585,   28711,   28835,   28955,      29073,   29188,   29299,
 245         29409,   29515,   29619,   29720,   29818,      29914,   30007,   30098,
 246         30186,   30272,   30356,   30437,   30516,      30593,   30668,   30740,
 247         30810,   30879,   30945,   31010,   31072,      31133,   31192,   31249,
 248
 249         31304,   31358,   31410,   31460,   31509,      31556,   31601,   31646,
 250         31688,   31730,   31770,   31808,   31846,      31882,   31917,   31950,
 251         31983,   32014,   32044,   32074,   32102,      32129,   32155,   32180,
 252         32205,   32228,   32251,   32273,   32294,      32314,   32333,   32352,
 253         32370,   32387,   32404,   32420,   32435,      32450,   32464,   32477,
 254         32490,   32503,   32515,   32526,   32537,      32548,   32558,   32568,
 255         32577,   32586,   32595,   32603,   32611,      32618,   32625,   32632,
 256         32639,   32645,   32651,   32657,   32662,      32667,   32672,   32677,
 257
 258         32682,   32686,   32690,   32694,   32698,      32702,   32705,   32708,
 259         32711,   32714,   32717,   32720,   32722,      32725,   32727,   32729,
 260         32731,   32733,   32735,   32737,   32739,      32740,   32742,   32743,
 261         32745,   32746,   32747,   32748,   32749,      32750,   32751,   32752,
 262         32753,   32754,   32755,   32756,   32757,      32757,   32758,   32758,
 263         32759,   32760,   32760,   32761,   32761,      32761,   32762,   32762,
 264         32763,   32763,   32763,   32764,   32764,      32764,   32764,   32765,
 265         32765,   32765,   32765,   32766,   32766,      32766,   32766,   32767,
 266 };
 267
 268
 269
 270 /*
 271  * Generate a random integer number of NORMAL distribution
 272  *
 273  * The table above is used to generate a pseudo-normal distribution,
 274  * in a manner which is much faster than calling a transcendental
 275  * function to calculate a true normal distribution.
 276  *
 277  * Basically, entry 64*N in the table above represents the number of
 278  * times out of 32767 that a random variable with normal distribution
 279  * will fall within N standard deviations of the mean.  That is, about
 280  * 68 percent of the time for N=1 and 95 percent of the time for N=2.
 281  *
 282  * The table above contains a "faked" final entry which allows us to
 283  * pretend that all values in a normal distribution are strictly less
 284  * than four standard deviations away from the mean.  This results in
 285  * "conservative" distribution of approximately 1/32768 values.
 286  *
 287  * Note that the binary search takes up to 16 quick iterations.
 288  */
 289 s16b randnor(int mean, int stand)
 290 {
 291         s16b tmp;
 292         s16b offset;
 293
 294         s16b low = 0;
 295         s16b high = RANDNOR_NUM;
 296         if (stand < 1) return (s16b)(mean);
 297
 298         /* Roll for probability */
 299         tmp = (s16b)randint0(32768);
 300
 301         /* Binary Search */
 302         while (low < high)
 303         {
 304                 int mid = (low + high) >> 1;
 305
 306                 /* Move right if forced */
 307                 if (randnor_table[mid] < tmp)
 308                 {
 309                         low = mid + 1;
 310                 }
 311
 312                 /* Move left otherwise */
 313                 else
 314                 {
 315                         high = (s16b)mid;
 316                 }
 317         }
 318
 319         /* Convert the index into an offset */
 320         offset = (long)stand * (long)low / RANDNOR_STD;
 321
 322         /* One half should be negative */
 323         if (randint0(100) < 50) return (mean - offset);
 324
 325         /* One half should be positive */
 326         return (mean + offset);
 327 }
 328
 329
 330
 331 /*
 332  * Generates damage for "2d6" style dice rolls
 333  */
 334 s16b damroll(DICE_NUMBER num, DICE_SID sides)
 335 {
 336         int i, sum = 0;
 337         for (i = 0; i < num; i++) sum += randint1(sides);
 338         return (s16b)(sum);
 339 }
 340
 341
 342 /*
 343  * Same as above, but always maximal
 344  */
 345 s16b maxroll(DICE_NUMBER num, DICE_SID sides)
 346 {
 347         return (num * sides);
 348 }
 349
 350
 351 /*
 352  * Given a numerator and a denominator, supply a properly rounded result,
 353  * using the RNG to smooth out remainders.  -LM-
 354  */
 355 s32b div_round(s32b n, s32b d)
 356 {
 357         s32b tmp;
 358
 359         /* Refuse to divide by zero */
 360         if (!d) return (n);
 361
 362         /* Division */
 363         tmp = n / d;
 364
 365         /* Rounding */
 366         if ((ABS(n) % ABS(d)) > randint0(ABS(d)))
 367         {
 368                 /* Increase the absolute value */
 369                 if (n * d > 0L) tmp += 1L;
 370                 else            tmp -= 1L;
 371         }
 372
 373         /* Return */
 374         return (tmp);
 375 }
 376
 377
 378
 379
 380 /*
 381  * Extract a "random" number from 0 to m-1, using the RNG.
 382  *
 383  * This function should be used when generating random numbers in
 384  * "external" program parts like the main-*.c files.  It preserves
 385  * the current RNG state to prevent influences on game-play.
 386  *
 387  * Could also use rand() from <stdlib.h> directly.
 388  */
 389 s32b Rand_external(s32b m)
 390 {
 391         static bool initialized = FALSE;
 392         static u32b Rand_state_external[4];
 393
 394         if (!initialized)
 395         {
 396                 /* Initialize with new seed */
 397                 u32b seed = (u32b)time(NULL);
 398                 Rand_seed(seed, Rand_state_external);
 399                 initialized = TRUE;
 400         }
 401
 402         return Rand_div_impl(m, Rand_state_external);
 403 }