1 |
48 |
robfinch |
/**********************************************************************
|
2 |
|
|
Freeciv - Copyright (C) 1996 - A Kjeldberg, L Gregersen, P Unold
|
3 |
|
|
This program is free software; you can redistribute it and/or modify
|
4 |
|
|
it under the terms of the GNU General Public License as published by
|
5 |
|
|
the Free Software Foundation; either version 2, or (at your option)
|
6 |
|
|
any later version.
|
7 |
|
|
|
8 |
|
|
This program is distributed in the hope that it will be useful,
|
9 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
10 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
11 |
|
|
GNU General Public License for more details.
|
12 |
|
|
***********************************************************************/
|
13 |
|
|
|
14 |
|
|
/*************************************************************************
|
15 |
|
|
The following random number generator can be found in _The Art of
|
16 |
|
|
Computer Programming Vol 2._ (2nd ed) by Donald E. Knuth. (C) 1998.
|
17 |
|
|
The algorithm is described in section 3.2.2 as Mitchell and Moore's
|
18 |
|
|
variant of a standard additive number generator. Note that the
|
19 |
|
|
the constants 55 and 24 are not random. Please become familiar with
|
20 |
|
|
this algorithm before you mess with it.
|
21 |
|
|
|
22 |
|
|
Since the additive number generator requires a table of numbers from
|
23 |
|
|
which to generate its random sequences, we must invent a way to
|
24 |
|
|
populate that table from a single seed value. I have chosen to do
|
25 |
|
|
this with a different PRNG, known as the "linear congruential method"
|
26 |
|
|
(also found in Knuth, Vol2). I must admit that my choices of constants
|
27 |
|
|
(3, 257, and MAX_UINT32) are probably not optimal, but they seem to
|
28 |
|
|
work well enough for our purposes.
|
29 |
|
|
|
30 |
|
|
Original author for this code: Cedric Tefft <cedric@earthling.net>
|
31 |
|
|
Modified to use rand_state struct by David Pfitzner <dwp@mso.anu.edu.au>
|
32 |
|
|
*************************************************************************/
|
33 |
|
|
|
34 |
|
|
#include "stdafx.h"
|
35 |
|
|
|
36 |
|
|
RANDOM_STATE RTFClasses::Random::rand_state;
|
37 |
|
|
|
38 |
|
|
/*************************************************************************
|
39 |
|
|
Returns a new random value from the sequence, in the interval 0 to
|
40 |
|
|
(size-1) inclusive, and updates global state for next call.
|
41 |
|
|
This means that if size <= 1 the function will always return 0.
|
42 |
|
|
|
43 |
|
|
Once we calculate new_rand below uniform (we hope) between 0 and
|
44 |
|
|
MAX_UINT32 inclusive, need to reduce to required range. Using
|
45 |
|
|
modulus is bad because generators like this are generally less
|
46 |
|
|
random for their low-significance bits, so this can give poor
|
47 |
|
|
results when 'size' is small. Instead want to divide the range
|
48 |
|
|
0..MAX_UINT32 into (size) blocks, each with (divisor) values, and
|
49 |
|
|
for any remainder, repeat the calculation of new_rand.
|
50 |
|
|
Then:
|
51 |
|
|
return_val = new_rand / divisor;
|
52 |
|
|
Will repeat for new_rand > max, where:
|
53 |
|
|
max = size * divisor - 1
|
54 |
|
|
Then max <= MAX_UINT32 implies
|
55 |
|
|
size * divisor <= (MAX_UINT32+1)
|
56 |
|
|
thus divisor <= (MAX_UINT32+1)/size
|
57 |
|
|
|
58 |
|
|
Need to calculate this divisor. Want divisor as large as possible
|
59 |
|
|
given above contraint, but it doesn't hurt us too much if it is a
|
60 |
|
|
bit smaller (just have to repeat more often). Calculation exactly
|
61 |
|
|
as above is complicated by fact that (MAX_UINT32+1) may not be
|
62 |
|
|
directly representable in type RANDOM_TYPE, so we do instead:
|
63 |
|
|
divisor = MAX_UINT32/size
|
64 |
|
|
*************************************************************************/
|
65 |
|
|
namespace RTFClasses
|
66 |
|
|
{
|
67 |
|
|
RANDOM_TYPE Random::rand(RANDOM_TYPE size)
|
68 |
|
|
{
|
69 |
|
|
RANDOM_TYPE new_rand, divisor, max;
|
70 |
|
|
int bailout = 0;
|
71 |
|
|
|
72 |
|
|
if (size > 1) {
|
73 |
|
|
divisor = MAX_UINT32 / size;
|
74 |
|
|
max = size * divisor - 1;
|
75 |
|
|
}
|
76 |
|
|
else {
|
77 |
|
|
/* size == 0 || size == 1 */
|
78 |
|
|
|
79 |
|
|
/*
|
80 |
|
|
* These assignments are only here to make the compiler
|
81 |
|
|
* happy. Since each usage is guarded with a if(size>1).
|
82 |
|
|
*/
|
83 |
|
|
max = MAX_UINT32;
|
84 |
|
|
divisor = 1;
|
85 |
|
|
}
|
86 |
|
|
|
87 |
|
|
do {
|
88 |
|
|
new_rand = (rand_state.v[rand_state.j]
|
89 |
|
|
+ rand_state.v[rand_state.k]) & MAX_UINT32;
|
90 |
|
|
|
91 |
|
|
rand_state.x = (rand_state.x +1) % 56;
|
92 |
|
|
rand_state.j = (rand_state.j +1) % 56;
|
93 |
|
|
rand_state.k = (rand_state.k +1) % 56;
|
94 |
|
|
rand_state.v[rand_state.x] = new_rand;
|
95 |
|
|
|
96 |
|
|
if (++bailout > 10000) {
|
97 |
|
|
// freelog(LOG_ERROR, "Bailout in myrand(%u)", size);
|
98 |
|
|
new_rand = 0;
|
99 |
|
|
break;
|
100 |
|
|
}
|
101 |
|
|
|
102 |
|
|
}
|
103 |
|
|
while (size > 1 && new_rand > max);
|
104 |
|
|
|
105 |
|
|
if (size > 1)
|
106 |
|
|
new_rand /= divisor;
|
107 |
|
|
else
|
108 |
|
|
new_rand = 0;
|
109 |
|
|
|
110 |
|
|
/* freelog(LOG_DEBUG, "rand(%u) = %u", size, new_rand); */
|
111 |
|
|
|
112 |
|
|
return new_rand;
|
113 |
|
|
}
|
114 |
|
|
|
115 |
|
|
/*************************************************************************
|
116 |
|
|
Initialize the generator; see comment at top of file.
|
117 |
|
|
*************************************************************************/
|
118 |
|
|
void Random::srand(RANDOM_TYPE seed)
|
119 |
|
|
{
|
120 |
|
|
int i;
|
121 |
|
|
|
122 |
|
|
rand_state.v[0]=(seed & MAX_UINT32);
|
123 |
|
|
|
124 |
|
|
for(i=1; i<56; i++)
|
125 |
|
|
rand_state.v[i] = (3 * rand_state.v[i-1] + 257) & MAX_UINT32;
|
126 |
|
|
|
127 |
|
|
rand_state.j = (55-55);
|
128 |
|
|
rand_state.k = (55-24);
|
129 |
|
|
rand_state.x = (55-0);
|
130 |
|
|
|
131 |
|
|
rand_state.is_init = true;
|
132 |
|
|
|
133 |
|
|
/* Heat it up a bit:
|
134 |
|
|
* Using modulus in myrand() this was important to pass
|
135 |
|
|
* test_random1(). Now using divisor in myrand() that particular
|
136 |
|
|
* test no longer indicates problems, but this seems a good idea
|
137 |
|
|
* anyway -- eg, other tests could well reveal other initial
|
138 |
|
|
* problems even using divisor.
|
139 |
|
|
*/
|
140 |
|
|
for (i=0; i<10000; i++)
|
141 |
|
|
(void) rand(MAX_UINT32);
|
142 |
|
|
}
|
143 |
|
|
|
144 |
|
|
/*************************************************************************
|
145 |
|
|
Test one aspect of randomness, using n numbers.
|
146 |
|
|
Reports results to LOG_NORMAL; with good randomness, behaviourchange
|
147 |
|
|
and behavioursame should be about the same size.
|
148 |
|
|
Tests current random state; saves and restores state, so can call
|
149 |
|
|
without interrupting current sequence.
|
150 |
|
|
*************************************************************************/
|
151 |
|
|
void Random::test(int n)
|
152 |
|
|
{
|
153 |
|
|
RANDOM_STATE saved_state;
|
154 |
|
|
int i, old_value = 0, new_value;
|
155 |
|
|
bool didchange, olddidchange = false;
|
156 |
|
|
int behaviourchange = 0, behavioursame = 0;
|
157 |
|
|
|
158 |
|
|
saved_state = getRandState();
|
159 |
|
|
/* mysrand(time(NULL)); */ /* use current state */
|
160 |
|
|
|
161 |
|
|
for (i = 0; i < n+2; i++) {
|
162 |
|
|
new_value = rand(2);
|
163 |
|
|
if (i > 0) { /* have old */
|
164 |
|
|
didchange = (new_value != old_value);
|
165 |
|
|
if (i > 1) { /* have olddidchange */
|
166 |
|
|
if (didchange != olddidchange)
|
167 |
|
|
behaviourchange++;
|
168 |
|
|
else
|
169 |
|
|
behavioursame++;
|
170 |
|
|
}
|
171 |
|
|
olddidchange = didchange;
|
172 |
|
|
}
|
173 |
|
|
old_value = new_value;
|
174 |
|
|
}
|
175 |
|
|
|
176 |
|
|
/* restore state: */
|
177 |
|
|
setRandState(saved_state);
|
178 |
|
|
};
|
179 |
|
|
|
180 |
|
|
}
|