1 /**
2 * \file z-rand.c
3 * \brief A Random Number Generator for Angband
4 *
5 * Copyright (c) 1997 Ben Harrison, Randy Hutson
6 *
7 * See below for copyright on the WELL random number generator.
8 *
9 * This work is free software; you can redistribute it and/or modify it
10 * under the terms of either:
11 *
12 * a) the GNU General Public License as published by the Free Software
13 * Foundation, version 2, or
14 *
15 * b) the "Angband licence":
16 * This software may be copied and distributed for educational, research,
17 * and not for profit purposes provided that this copyright and statement
18 * are included in all such copies. Other copyrights may also apply.
19 */
20 #include "z-rand.h"
21
22 /**
23 * This file provides a pseudo-random number generator.
24 *
25 * This code provides both a "quick" random number generator (4 bytes of
26 * state), and a "complex" random number generator (128 + 4 bytes of state).
27 *
28 * The complex RNG (used for most game entropy) is provided by the WELL102a
29 * algorithm, used with permission. See below for copyright information
30 * about the WELL implementation.
31 *
32 * To use of the "simple" RNG, activate it via "Rand_quick = true" and
33 * "Rand_value = seed". After that it will be automatically used instead of
34 * the "complex" RNG. When you are done, you can de-activate it via
35 * "Rand_quick = false". You can also choose a new seed.
36 */
37
38 /* begin WELL RNG
39 * *************************************************************************
40 * Copyright: Francois Panneton and Pierre L'Ecuyer, University of Montreal
41 * Makoto Matsumoto, Hiroshima University
42 * *************************************************************************
43 * Code was modified slightly by Erik Osheim to work on unsigned integers.
44 */
45 #define M1 3
46 #define M2 24
47 #define M3 10
48
49 #define MAT0POS(t, v) (v ^ (v >> t))
50 #define MAT0NEG(t, v) (v ^ (v << (-(t))))
51 #define Identity(v) (v)
52
53 u32b state_i = 0;
54 u32b STATE[RAND_DEG] = {0, 0, 0, 0, 0, 0, 0, 0,
55 0, 0, 0, 0, 0, 0, 0, 0,
56 0, 0, 0, 0, 0, 0, 0, 0,
57 0, 0, 0, 0, 0, 0, 0, 0};
58 u32b z0, z1, z2;
59
60 #define V0 STATE[state_i]
61 #define VM1 STATE[(state_i + M1) & 0x0000001fU]
62 #define VM2 STATE[(state_i + M2) & 0x0000001fU]
63 #define VM3 STATE[(state_i + M3) & 0x0000001fU]
64 #define VRm1 STATE[(state_i + 31) & 0x0000001fU]
65 #define newV0 STATE[(state_i + 31) & 0x0000001fU]
66 #define newV1 STATE[state_i]
67
WELLRNG1024a(void)68 static u32b WELLRNG1024a (void){
69 z0 = VRm1;
70 z1 = Identity(V0) ^ MAT0POS (8, VM1);
71 z2 = MAT0NEG (-19, VM2) ^ MAT0NEG(-14,VM3);
72 newV1 = z1 ^ z2;
73 newV0 = MAT0NEG (-11,z0) ^ MAT0NEG(-7,z1) ^ MAT0NEG(-13,z2);
74 state_i = (state_i + 31) & 0x0000001fU;
75 return STATE[state_i];
76 }
77 /* end WELL RNG */
78
79 /**
80 * Simple RNG, implemented with a linear congruent algorithm.
81 */
82 #define LCRNG(X) ((X) * 1103515245 + 12345)
83
84
85 /**
86 * Whether to use the simple RNG or not.
87 */
88 bool Rand_quick = true;
89
90 /**
91 * The current "seed" of the simple RNG.
92 */
93 u32b Rand_value;
94
95 static bool rand_fixed = false;
96 static u32b rand_fixval = 0;
97
98 /**
99 * Initialize the complex RNG using a new seed.
100 */
Rand_state_init(u32b seed)101 void Rand_state_init(u32b seed)
102 {
103 int i, j;
104
105 /* Seed the table */
106 STATE[0] = seed;
107
108 /* Propagate the seed */
109 for (i = 1; i < RAND_DEG; i++)
110 STATE[i] = LCRNG(STATE[i - 1]);
111
112 /* Cycle the table ten times per degree */
113 for (i = 0; i < RAND_DEG * 10; i++) {
114 /* Acquire the next index */
115 j = (state_i + 1) % RAND_DEG;
116
117 /* Update the table, extract an entry */
118 STATE[j] += STATE[state_i];
119
120 /* Advance the index */
121 state_i = j;
122 }
123 }
124
125 /**
126 * Initialise the RNG
127 */
Rand_init(void)128 void Rand_init(void)
129 {
130 /* Init RNG */
131 if (Rand_quick) {
132 u32b seed;
133
134 /* Basic seed */
135 seed = (u32b)(time(NULL));
136
137 #ifdef UNIX
138
139 /* Mutate the seed on Unix machines */
140 seed = ((seed >> 3) * (getpid() << 1));
141
142 #endif
143
144 /* Use the complex RNG */
145 Rand_quick = false;
146
147 /* Seed the "complex" RNG */
148 Rand_state_init(seed);
149 }
150 }
151
152
153 /**
154 * Extract a "random" number from 0 to m - 1, via division.
155 *
156 * This method selects "random" 28-bit numbers, and then uses division to drop
157 * those numbers into "m" different partitions, plus a small non-partition to
158 * reduce bias, taking as the final value the first "good" partition that a
159 * number falls into.
160 *
161 * This method has no bias, and is much less affected by patterns in the "low"
162 * bits of the underlying RNG's. However, it is potentially non-terminating.
163 */
Rand_div(u32b m)164 u32b Rand_div(u32b m)
165 {
166 u32b n, r = 0;
167
168 /* Division by zero will result if m is larger than 0x10000000 */
169 assert(m <= 0x10000000);
170
171 /* Hack -- simple case */
172 if (m <= 1) return (0);
173
174 if (rand_fixed)
175 return (rand_fixval * 1000 * (m - 1)) / (100 * 1000);
176
177 /* Partition size */
178 n = (0x10000000 / m);
179
180 if (Rand_quick) {
181 /* Use a simple RNG */
182 /* Wait for it */
183 while (1) {
184 /* Cycle the generator */
185 r = (Rand_value = LCRNG(Rand_value));
186
187 /* Mutate a 28-bit "random" number */
188 r = ((r >> 4) & 0x0FFFFFFF) / n;
189
190 /* Done */
191 if (r < m) break;
192 }
193 } else {
194 /* Use a complex RNG */
195 while (1) {
196 /* Get the next pseudorandom number */
197 r = WELLRNG1024a();
198
199 /* Mutate a 28-bit "random" number */
200 r = ((r >> 4) & 0x0FFFFFFF) / n;
201
202 /* Done */
203 if (r < m) break;
204 }
205 }
206
207 /* Use the value */
208 return (r);
209 }
210
211
212 /**
213 * The number of entries in the "Rand_normal_table"
214 */
215 #define RANDNOR_NUM 256
216
217 /**
218 * The standard deviation of the "Rand_normal_table"
219 */
220 #define RANDNOR_STD 64
221
222 /**
223 * The normal distribution table for the "Rand_normal()" function (below)
224 */
225 static s16b Rand_normal_table[RANDNOR_NUM] = {
226 206, 613, 1022, 1430, 1838, 2245, 2652, 3058,
227 3463, 3867, 4271, 4673, 5075, 5475, 5874, 6271,
228 6667, 7061, 7454, 7845, 8234, 8621, 9006, 9389,
229 9770, 10148, 10524, 10898, 11269, 11638, 12004, 12367,
230 12727, 13085, 13440, 13792, 14140, 14486, 14828, 15168,
231 15504, 15836, 16166, 16492, 16814, 17133, 17449, 17761,
232 18069, 18374, 18675, 18972, 19266, 19556, 19842, 20124,
233 20403, 20678, 20949, 21216, 21479, 21738, 21994, 22245,
234
235 22493, 22737, 22977, 23213, 23446, 23674, 23899, 24120,
236 24336, 24550, 24759, 24965, 25166, 25365, 25559, 25750,
237 25937, 26120, 26300, 26476, 26649, 26818, 26983, 27146,
238 27304, 27460, 27612, 27760, 27906, 28048, 28187, 28323,
239 28455, 28585, 28711, 28835, 28955, 29073, 29188, 29299,
240 29409, 29515, 29619, 29720, 29818, 29914, 30007, 30098,
241 30186, 30272, 30356, 30437, 30516, 30593, 30668, 30740,
242 30810, 30879, 30945, 31010, 31072, 31133, 31192, 31249,
243
244 31304, 31358, 31410, 31460, 31509, 31556, 31601, 31646,
245 31688, 31730, 31770, 31808, 31846, 31882, 31917, 31950,
246 31983, 32014, 32044, 32074, 32102, 32129, 32155, 32180,
247 32205, 32228, 32251, 32273, 32294, 32314, 32333, 32352,
248 32370, 32387, 32404, 32420, 32435, 32450, 32464, 32477,
249 32490, 32503, 32515, 32526, 32537, 32548, 32558, 32568,
250 32577, 32586, 32595, 32603, 32611, 32618, 32625, 32632,
251 32639, 32645, 32651, 32657, 32662, 32667, 32672, 32677,
252
253 32682, 32686, 32690, 32694, 32698, 32702, 32705, 32708,
254 32711, 32714, 32717, 32720, 32722, 32725, 32727, 32729,
255 32731, 32733, 32735, 32737, 32739, 32740, 32742, 32743,
256 32745, 32746, 32747, 32748, 32749, 32750, 32751, 32752,
257 32753, 32754, 32755, 32756, 32757, 32757, 32758, 32758,
258 32759, 32760, 32760, 32761, 32761, 32761, 32762, 32762,
259 32763, 32763, 32763, 32764, 32764, 32764, 32764, 32765,
260 32765, 32765, 32765, 32766, 32766, 32766, 32766, 32767,
261 };
262
263
264 /**
265 * Generate a random integer number of NORMAL distribution
266 *
267 * The table above is used to generate a psuedo-normal distribution, in a
268 * manner which is much faster than calling a transcendental function to
269 * calculate a true normal distribution.
270 *
271 * Basically, entry 64 * N in the table above represents the number of times
272 * out of 32767 that a random variable with normal distribution will fall
273 * within N standard deviations of the mean. That is, about 68 percent of the
274 * time for N=1 and 95 percent of the time for N=2.
275 *
276 * The table above contains a "faked" final entry which allows us to pretend
277 * that all values in a normal distribution are strictly less than four
278 * standard deviations away from the mean. This results in "conservative"
279 * distribution of approximately 1/32768 values.
280 *
281 * Note that the binary search takes up to 16 quick iterations.
282 */
Rand_normal(int mean,int stand)283 s16b Rand_normal(int mean, int stand)
284 {
285 s16b tmp, offset;
286
287 s16b low = 0;
288 s16b high = RANDNOR_NUM;
289
290 /* Paranoia */
291 if (stand < 1) return (mean);
292
293 /* Roll for probability */
294 tmp = (s16b)randint0(32768);
295
296 /* Binary Search */
297 while (low < high) {
298 int mid = (low + high) >> 1;
299
300 /* Move right if forced */
301 if (Rand_normal_table[mid] < tmp) {
302 low = mid + 1;
303 } else {
304 high = mid;
305 }
306 }
307
308 /* Convert the index into an offset */
309 offset = (s16b)((long)stand * (long)low / RANDNOR_STD);
310
311 /* One half should be negative */
312 if (one_in_(2)) return (mean - offset);
313
314 /* One half should be positive */
315 return (mean + offset);
316 }
317
318
319 /**
320 * Choose an integer from a distribution where we know the mean and approximate
321 * upper and lower bounds.
322 *
323 * We divide the imagined distribution into two halves, above and below the
324 * mean, and then treat the bounds as if they are the given number of
325 * standard deviations from the mean in the appropriate direction. Note that
326 * `stand_u` and `stand_l` are 10 times the number of standart deviations we
327 * are asking for.
328 * The function chooses an integer from a normal distribution, and then scales
329 * it to fit the target distribution.
330 */
Rand_sample(int mean,int upper,int lower,int stand_u,int stand_l)331 int Rand_sample(int mean, int upper, int lower, int stand_u, int stand_l)
332 {
333 int pick = Rand_normal(0, 1000);
334
335 /* Scale to fit */
336 if (pick > 0) {
337 /* Positive pick, scale up */
338 pick *= (upper - mean);
339 pick /= (100 * stand_u);
340 } else if (pick < 0) {
341 /* Negative pick, scale down */
342 pick *= (mean - lower);
343 pick /= (100 * stand_l);
344 }
345
346 return mean + pick;
347 }
348
349 /**
350 * Generates damage for "2d6" style dice rolls
351 */
damroll(int num,int sides)352 int damroll(int num, int sides)
353 {
354 int i;
355 int sum = 0;
356
357 if (sides <= 0) return 0;
358
359 for (i = 0; i < num; i++)
360 sum += randint1(sides);
361 return sum;
362 }
363
364
365
366 /**
367 * Calculation helper function for damroll
368 */
damcalc(int num,int sides,aspect dam_aspect)369 int damcalc(int num, int sides, aspect dam_aspect)
370 {
371 switch (dam_aspect) {
372 case MAXIMISE:
373 case EXTREMIFY: return num * sides;
374 case RANDOMISE: return damroll(num, sides);
375 case MINIMISE: return num;
376 case AVERAGE: return num * (sides + 1) / 2;
377 }
378
379 assert(0 && "Should never reach here");
380 return 0;
381 }
382
383
384 /**
385 * Generates a random signed long integer X where `A` <= X <= `B`.
386 * The integer X falls along a uniform distribution.
387 *
388 * Note that "rand_range(0, N-1)" == "randint0(N)".
389 */
rand_range(int A,int B)390 int rand_range(int A, int B)
391 {
392 if (A == B) return A;
393 assert(A < B);
394
395 return A + (s32b)Rand_div(1 + B - A);
396 }
397
398
399 /**
400 * Perform division, possibly rounding up or down depending on the size of the
401 * remainder and chance.
402 */
simulate_division(int dividend,int divisor)403 static int simulate_division(int dividend, int divisor)
404 {
405 int quotient = dividend / divisor;
406 int remainder = dividend % divisor;
407 if (randint0(divisor) < remainder) quotient++;
408 return quotient;
409 }
410
411
412 /**
413 * Help determine an "enchantment bonus" for an object.
414 *
415 * To avoid floating point but still provide a smooth distribution of bonuses,
416 * we simply round the results of division in such a way as to "average" the
417 * correct floating point value.
418 *
419 * This function has been changed. It uses "Rand_normal()" to choose values
420 * from a normal distribution, whose mean moves from zero towards the max as
421 * the level increases, and whose standard deviation is equal to 1/4 of the
422 * max, and whose values are forced to lie between zero and the max, inclusive.
423 *
424 * Since the "level" rarely passes 100 before Morgoth is dead, it is very
425 * rare to get the "full" enchantment on an object, even a deep levels.
426 *
427 * It is always possible (albeit unlikely) to get the "full" enchantment.
428 *
429 * A sample distribution of values from "m_bonus(10, N)" is shown below:
430 *
431 * N 0 1 2 3 4 5 6 7 8 9 10
432 * --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
433 * 0 66.37 13.01 9.73 5.47 2.89 1.31 0.72 0.26 0.12 0.09 0.03
434 * 8 46.85 24.66 12.13 8.13 4.20 2.30 1.05 0.36 0.19 0.08 0.05
435 * 16 30.12 27.62 18.52 10.52 6.34 3.52 1.95 0.90 0.31 0.15 0.05
436 * 24 22.44 15.62 30.14 12.92 8.55 5.30 2.39 1.63 0.62 0.28 0.11
437 * 32 16.23 11.43 23.01 22.31 11.19 7.18 4.46 2.13 1.20 0.45 0.41
438 * 40 10.76 8.91 12.80 29.51 16.00 9.69 5.90 3.43 1.47 0.88 0.65
439 * 48 7.28 6.81 10.51 18.27 27.57 11.76 7.85 4.99 2.80 1.22 0.94
440 * 56 4.41 4.73 8.52 11.96 24.94 19.78 11.06 7.18 3.68 1.96 1.78
441 * 64 2.81 3.07 5.65 9.17 13.01 31.57 13.70 9.30 6.04 3.04 2.64
442 * 72 1.87 1.99 3.68 7.15 10.56 20.24 25.78 12.17 7.52 4.42 4.62
443 * 80 1.02 1.23 2.78 4.75 8.37 12.04 27.61 18.07 10.28 6.52 7.33
444 * 88 0.70 0.57 1.56 3.12 6.34 10.06 15.76 30.46 12.58 8.47 10.38
445 * 96 0.27 0.60 1.25 2.28 4.30 7.60 10.77 22.52 22.51 11.37 16.53
446 * 104 0.22 0.42 0.77 1.36 2.62 5.33 8.93 13.05 29.54 15.23 22.53
447 * 112 0.15 0.20 0.56 0.87 2.00 3.83 6.86 10.06 17.89 27.31 30.27
448 * 120 0.03 0.11 0.31 0.46 1.31 2.48 4.60 7.78 11.67 25.53 45.72
449 * 128 0.02 0.01 0.13 0.33 0.83 1.41 3.24 6.17 9.57 14.22 64.07
450 */
m_bonus(int max,int level)451 s16b m_bonus(int max, int level)
452 {
453 int bonus, stand, value;
454
455 /* Make sure level is reasonable */
456 if (level >= MAX_RAND_DEPTH) level = MAX_RAND_DEPTH - 1;
457
458 /* The bonus approaches max as level approaches MAX_RAND_DEPTH */
459 bonus = simulate_division(max * level, MAX_RAND_DEPTH);
460
461 /* The standard deviation is 1/4 of the max */
462 stand = simulate_division(max, 4);
463
464 /* Choose a value */
465 value = Rand_normal(bonus, stand);
466
467 /* Return, enforcing the min and max values */
468 if (value < 0)
469 return 0;
470 else if (value > max)
471 return max;
472 else
473 return value;
474 }
475
476
477 /**
478 * Calculation helper function for m_bonus
479 */
m_bonus_calc(int max,int level,aspect bonus_aspect)480 s16b m_bonus_calc(int max, int level, aspect bonus_aspect)
481 {
482 switch (bonus_aspect) {
483 case EXTREMIFY:
484 case MAXIMISE: return max;
485 case RANDOMISE: return m_bonus(max, level);
486 case MINIMISE: return 0;
487 case AVERAGE: return max * level / MAX_RAND_DEPTH;
488 }
489
490 assert(0 && "Should never reach here");
491 return 0;
492 }
493
494
495 /**
496 * Calculation helper function for random_value structs
497 */
randcalc(random_value v,int level,aspect rand_aspect)498 int randcalc(random_value v, int level, aspect rand_aspect)
499 {
500 if (rand_aspect == EXTREMIFY) {
501 int min = randcalc(v, level, MINIMISE);
502 int max = randcalc(v, level, MAXIMISE);
503 return abs(min) > abs(max) ? min : max;
504
505 } else {
506 int dmg = damcalc(v.dice, v.sides, rand_aspect);
507 int bonus = m_bonus_calc(v.m_bonus, level, rand_aspect);
508 return v.base + dmg + bonus;
509 }
510 }
511
512
513 /**
514 * Test to see if a value is within a random_value's range
515 */
randcalc_valid(random_value v,int test)516 bool randcalc_valid(random_value v, int test)
517 {
518 if (test < randcalc(v, 0, MINIMISE))
519 return false;
520 else if (test > randcalc(v, 0, MAXIMISE))
521 return false;
522 else
523 return true;
524 }
525
526 /**
527 * Test to see if a random_value actually varies
528 */
randcalc_varies(random_value v)529 bool randcalc_varies(random_value v)
530 {
531 return randcalc(v, 0, MINIMISE) != randcalc(v, 0, MAXIMISE);
532 }
533
rand_fix(u32b val)534 void rand_fix(u32b val)
535 {
536 rand_fixed = true;
537 rand_fixval = val;
538 }
539
540 int getpid(void);
541
542 /**
543 * Another simple RNG that does not use any of the above state
544 * (so can be used without disturbing the game's RNG state)
545 */
Rand_simple(u32b m)546 u32b Rand_simple(u32b m)
547 {
548 static time_t seed;
549 time_t v;
550 v = time(NULL);
551 seed = LCRNG(seed % m) + ((v << 16) ^ v ^ getpid());
552 return (seed % m);
553 }
554