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 typedef struct {
  74         u32b dw[2];
  75 } u64b;
  76
  77 static u64b u64b_xor(u64b a, u64b b)
  78 {
  79         u64b result;
  80
  81         result.dw[0] = a.dw[0] ^ b.dw[0];
  82         result.dw[1] = a.dw[1] ^ b.dw[1];
  83
  84         return result;
  85 }
  86
  87 static u64b u64b_shiftl(u64b x, int k)
  88 {
  89         u64b result;
  90
  91         if (k < 32) {
  92                 result.dw[1] = (x.dw[1] << k) | (x.dw[0] >> (32 - k));
  93                 result.dw[0] = (x.dw[0] << k);
  94         }
  95         else {
  96                 result.dw[1] = (x.dw[0] << (k - 32));
  97                 result.dw[0] = 0;
  98         }
  99
 100         return result;
 101 }
 102
 103 static u64b u64b_rotl(u64b x, int k)
 104 {
 105         u64b result;
 106
 107         if (k < 32) {
 108                 result.dw[0] = (x.dw[0] << k) | (x.dw[1] >> (32 - k));
 109                 result.dw[1] = (x.dw[1] << k) | (x.dw[0] >> (32 - k));
 110         }
 111         else {
 112                 result.dw[0] = (x.dw[0] >> (64 - k)) | (x.dw[1] << (k - 32));
 113                 result.dw[1] = (x.dw[1] >> (64 - k)) | (x.dw[0] << (k - 32));
 114         }
 115
 116         return result;
 117 }
 118
 119 static u64b u64b_add(u64b a, u64b b)
 120 {
 121         u64b result;
 122
 123         result.dw[0] = a.dw[0] + b.dw[0];
 124         result.dw[1] = a.dw[1] + b.dw[1];
 125
 126         if (result.dw[0] < a.dw[0])
 127                 result.dw[1] ++;
 128
 129         return result;
 130 }
 131
 132 /*
 133  * Initialize Xorshift Algorithm state
 134  */
 135 static void Rand_Xorshift_seed(u32b seed, u32b* state)
 136 {
 137         int i;
 138
 139         /* Initialize Xorshift Algorithm RNG */
 140         for (i = 1; i <= 4; ++ i) {
 141                 seed = 1812433253UL * (seed ^ (seed >> 30)) + i;
 142                 state[i-1] = seed;
 143         }
 144 }
 145
 146 #if 0
 147 /*
 148  * Xorshift Algorithm
 149  */
 150 static u32b Rand_Xorshift(u32b* state)
 151 {
 152         u32b t = state[0] ^ (state[0] << 11);
 153
 154         state[0] = state[1];
 155         state[1] = state[2];
 156         state[2] = state[3];
 157
 158         state[3] = (state[3] ^ (state[3] >> 19)) ^ (t ^ (t >> 8));
 159
 160         return state[3];
 161 }
 162 #endif
 163
 164 /*
 165  * Xoroshiro128+ Algorithm
 166  */
 167 static u32b Rand_Xoroshiro128plus(u32b* state)
 168 {
 169         const u64b s0 = *((u64b*)state);
 170         u64b s1 = *((u64b*)state + 1);
 171         const u64b result = u64b_add(s0, s1);
 172
 173         s1 = u64b_xor(s0, s1);
 174         *((u64b*)state) = u64b_xor(u64b_xor(u64b_rotl(s0, 55), s1), u64b_shiftl(s1, 14));
 175         *((u64b*)state + 1) = u64b_rotl(s1, 36);
 176
 177         return result.dw[0];
 178 }
 179
 180 static const u32b Rand_Xorshift_max = 0xFFFFFFFF;
 181
 182 /*
 183  * Initialize the RNG using a new seed
 184  */
 185 void Rand_state_set(u32b seed)
 186 {
 187         Rand_Xorshift_seed(seed, Rand_state);
 188 }
 189
 190 void Rand_state_init(void)
 191 {
 192 #ifdef RNG_DEVICE
 193
 194         FILE *fp = fopen(RNG_DEVICE, "r");
 195
 196         do {
 197                 fread(Rand_state, sizeof(Rand_state[0]), 4, fp);
 198         } while ((Rand_state[0] | Rand_state[1] | Rand_state[2] | Rand_state[3]) == 0);
 199
 200         fclose(fp);
 201
 202 #elif defined(WINDOWS)
 203
 204         HCRYPTPROV hProvider;
 205
 206         CryptAcquireContext(&hProvider, NULL, NULL, PROV_RSA_FULL, 0);
 207
 208         do {
 209                 CryptGenRandom(hProvider, sizeof(Rand_state[0]) * 4, (BYTE*)Rand_state);
 210         } while ((Rand_state[0] | Rand_state[1] | Rand_state[2] | Rand_state[3]) == 0);
 211
 212         CryptReleaseContext(hProvider, 0);
 213
 214 #else
 215
 216         /* Basic seed */
 217         u32b seed = (time(NULL));
 218 #ifdef SET_UID
 219         /* Mutate the seed on Unix machines */
 220         seed = ((seed >> 3) * (getpid() << 1));
 221 #endif
 222         /* Seed the RNG */
 223         Rand_state_set(seed);
 224
 225 #endif
 226 }
 227
 228 /*
 229  * Backup the RNG state
 230  */
 231 void Rand_state_backup(u32b* backup_state)
 232 {
 233         int i;
 234
 235         for (i = 0; i < 4; ++ i) {
 236                 backup_state[i] = Rand_state[i];
 237         }
 238 }
 239
 240 /*
 241  * Restore the RNG state
 242  */
 243 void Rand_state_restore(u32b* backup_state)
 244 {
 245         int i;
 246
 247         for (i = 0; i < 4; ++ i) {
 248                 Rand_state[i] = backup_state[i];
 249         }
 250 }
 251
 252
 253 /*
 254  * Extract a "random" number from 0 to m-1, via "division"
 255  */
 256 static s32b Rand_div_impl(s32b m, u32b* state)
 257 {
 258         u32b scaling;
 259         u32b past;
 260         u32b ret;
 261
 262         /* Hack -- simple case */
 263         if (m <= 1) return (0);
 264
 265         scaling = Rand_Xorshift_max / m;
 266         past = scaling * m;
 267
 268         do {
 269                 ret = Rand_Xoroshiro128plus(state);
 270         } while (ret >= past);
 271
 272         return ret / scaling;
 273 }
 274
 275 s32b Rand_div(s32b m)
 276 {
 277         return Rand_div_impl(m, Rand_state);
 278 }
 279
 280
 281
 282
 283 /*
 284  * The number of entries in the "randnor_table"
 285  */
 286 #define RANDNOR_NUM     256
 287
 288 /*
 289  * The standard deviation of the "randnor_table"
 290  */
 291 #define RANDNOR_STD     64
 292
 293 /*
 294  * The normal distribution table for the "randnor()" function (below)
 295  */
 296 static s16b randnor_table[RANDNOR_NUM] =
 297 {
 298         206,     613,    1022,    1430,         1838,    2245,    2652,    3058,
 299         3463,    3867,    4271,    4673,        5075,    5475,    5874,    6271,
 300         6667,    7061,    7454,    7845,        8234,    8621,    9006,    9389,
 301         9770,   10148,   10524,   10898,   11269,       11638,   12004,   12367,
 302         12727,   13085,   13440,   13792,   14140,      14486,   14828,   15168,
 303         15504,   15836,   16166,   16492,   16814,      17133,   17449,   17761,
 304         18069,   18374,   18675,   18972,   19266,      19556,   19842,   20124,
 305         20403,   20678,   20949,   21216,   21479,      21738,   21994,   22245,
 306
 307         22493,   22737,   22977,   23213,   23446,      23674,   23899,   24120,
 308         24336,   24550,   24759,   24965,   25166,      25365,   25559,   25750,
 309         25937,   26120,   26300,   26476,   26649,      26818,   26983,   27146,
 310         27304,   27460,   27612,   27760,   27906,      28048,   28187,   28323,
 311         28455,   28585,   28711,   28835,   28955,      29073,   29188,   29299,
 312         29409,   29515,   29619,   29720,   29818,      29914,   30007,   30098,
 313         30186,   30272,   30356,   30437,   30516,      30593,   30668,   30740,
 314         30810,   30879,   30945,   31010,   31072,      31133,   31192,   31249,
 315
 316         31304,   31358,   31410,   31460,   31509,      31556,   31601,   31646,
 317         31688,   31730,   31770,   31808,   31846,      31882,   31917,   31950,
 318         31983,   32014,   32044,   32074,   32102,      32129,   32155,   32180,
 319         32205,   32228,   32251,   32273,   32294,      32314,   32333,   32352,
 320         32370,   32387,   32404,   32420,   32435,      32450,   32464,   32477,
 321         32490,   32503,   32515,   32526,   32537,      32548,   32558,   32568,
 322         32577,   32586,   32595,   32603,   32611,      32618,   32625,   32632,
 323         32639,   32645,   32651,   32657,   32662,      32667,   32672,   32677,
 324
 325         32682,   32686,   32690,   32694,   32698,      32702,   32705,   32708,
 326         32711,   32714,   32717,   32720,   32722,      32725,   32727,   32729,
 327         32731,   32733,   32735,   32737,   32739,      32740,   32742,   32743,
 328         32745,   32746,   32747,   32748,   32749,      32750,   32751,   32752,
 329         32753,   32754,   32755,   32756,   32757,      32757,   32758,   32758,
 330         32759,   32760,   32760,   32761,   32761,      32761,   32762,   32762,
 331         32763,   32763,   32763,   32764,   32764,      32764,   32764,   32765,
 332         32765,   32765,   32765,   32766,   32766,      32766,   32766,   32767,
 333 };
 334
 335
 336
 337 /*
 338  * Generate a random integer number of NORMAL distribution
 339  *
 340  * The table above is used to generate a pseudo-normal distribution,
 341  * in a manner which is much faster than calling a transcendental
 342  * function to calculate a true normal distribution.
 343  *
 344  * Basically, entry 64*N in the table above represents the number of
 345  * times out of 32767 that a random variable with normal distribution
 346  * will fall within N standard deviations of the mean.  That is, about
 347  * 68 percent of the time for N=1 and 95 percent of the time for N=2.
 348  *
 349  * The table above contains a "faked" final entry which allows us to
 350  * pretend that all values in a normal distribution are strictly less
 351  * than four standard deviations away from the mean.  This results in
 352  * "conservative" distribution of approximately 1/32768 values.
 353  *
 354  * Note that the binary search takes up to 16 quick iterations.
 355  */
 356 s16b randnor(int mean, int stand)
 357 {
 358         s16b tmp;
 359         s16b offset;
 360
 361         s16b low = 0;
 362         s16b high = RANDNOR_NUM;
 363
 364         /* Paranoia */
 365         if (stand < 1) return (s16b)(mean);
 366
 367         /* Roll for probability */
 368         tmp = (s16b)randint0(32768);
 369
 370         /* Binary Search */
 371         while (low < high)
 372         {
 373                 int mid = (low + high) >> 1;
 374
 375                 /* Move right if forced */
 376                 if (randnor_table[mid] < tmp)
 377                 {
 378                         low = mid + 1;
 379                 }
 380
 381                 /* Move left otherwise */
 382                 else
 383                 {
 384                         high = (s16b)mid;
 385                 }
 386         }
 387
 388         /* Convert the index into an offset */
 389         offset = (long)stand * (long)low / RANDNOR_STD;
 390
 391         /* One half should be negative */
 392         if (randint0(100) < 50) return (mean - offset);
 393
 394         /* One half should be positive */
 395         return (mean + offset);
 396 }
 397
 398
 399
 400 /*
 401  * Generates damage for "2d6" style dice rolls
 402  */
 403 s16b damroll(int num, int sides)
 404 {
 405         int i, sum = 0;
 406         for (i = 0; i < num; i++) sum += randint1(sides);
 407         return (s16b)(sum);
 408 }
 409
 410
 411 /*
 412  * Same as above, but always maximal
 413  */
 414 s16b maxroll(int num, int sides)
 415 {
 416         return (num * sides);
 417 }
 418
 419
 420 /*
 421  * Given a numerator and a denominator, supply a properly rounded result,
 422  * using the RNG to smooth out remainders.  -LM-
 423  */
 424 s32b div_round(s32b n, s32b d)
 425 {
 426         s32b tmp;
 427
 428         /* Refuse to divide by zero */
 429         if (!d) return (n);
 430
 431         /* Division */
 432         tmp = n / d;
 433
 434         /* Rounding */
 435         if ((ABS(n) % ABS(d)) > randint0(ABS(d)))
 436         {
 437                 /* Increase the absolute value */
 438                 if (n * d > 0L) tmp += 1L;
 439                 else            tmp -= 1L;
 440         }
 441
 442         /* Return */
 443         return (tmp);
 444 }
 445
 446
 447
 448
 449 /*
 450  * Extract a "random" number from 0 to m-1, using the RNG.
 451  *
 452  * This function should be used when generating random numbers in
 453  * "external" program parts like the main-*.c files.  It preserves
 454  * the current RNG state to prevent influences on game-play.
 455  *
 456  * Could also use rand() from <stdlib.h> directly. XXX XXX XXX
 457  */
 458 s32b Rand_external(s32b m)
 459 {
 460         static bool initialized = FALSE;
 461         static u32b Rand_state_external[4];
 462
 463         if (!initialized)
 464         {
 465                 /* Initialize with new seed */
 466                 u32b seed = (u32b)time(NULL);
 467                 Rand_Xorshift_seed(seed, Rand_state_external);
 468                 initialized = TRUE;
 469         }
 470
 471         return Rand_div_impl(m, Rand_state_external);
 472 }