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