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