1 #include <assert.h>
2 #include <math.h>
3 #include "simrandom.h"
4 #include "../simsys.h"
5
6 /* This is the mersenne random generator: More random and faster! */
7
8 /* Period parameters */
9 #define MERSENNE_TWISTER_N 624
10 #define M 397
11 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
12 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
13 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
14
15 static unsigned long mersenne_twister[MERSENNE_TWISTER_N]; // the array for the state vector
16 static int mersenne_twister_index = MERSENNE_TWISTER_N + 1; // mersenne_twister_index==N+1 means mersenne_twister[N] is not initialized
17
18 static uint8 random_origin = 0;
19
20
21 /* initializes mersenne_twister[N] with a seed */
init_genrand(uint32 s)22 static void init_genrand(uint32 s)
23 {
24 mersenne_twister[0]= s & 0xffffffffUL;
25 for (mersenne_twister_index=1; mersenne_twister_index<MERSENNE_TWISTER_N; mersenne_twister_index++) {
26 mersenne_twister[mersenne_twister_index] = (1812433253UL * (mersenne_twister[mersenne_twister_index-1] ^ (mersenne_twister[mersenne_twister_index-1] >> 30)) + mersenne_twister_index);
27 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
28 /* In the previous versions, MSBs of the seed affect */
29 /* only MSBs of the array mersenne_twister[]. */
30 /* 2002/01/09 modified by Makoto Matsumoto */
31 mersenne_twister[mersenne_twister_index] &= 0xffffffffUL;
32 /* for >32 bit machines */
33 }
34 }
35
36
37 /* generate N words at one time */
MTgenerate()38 static void MTgenerate()
39 {
40 static uint32 mag01[2]={0x0UL, MATRIX_A};
41 uint32 y;
42 int kk;
43
44 if (mersenne_twister_index == MERSENNE_TWISTER_N+1) /* if init_genrand() has not been called, */
45 init_genrand(5489UL); /* a default initial seed is used */
46
47 for (kk=0;kk<MERSENNE_TWISTER_N-M;kk++) {
48 y = (mersenne_twister[kk]&UPPER_MASK)|(mersenne_twister[kk+1]&LOWER_MASK);
49 mersenne_twister[kk] = mersenne_twister[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
50 }
51 for (;kk<MERSENNE_TWISTER_N-1;kk++) {
52 y = (mersenne_twister[kk]&UPPER_MASK)|(mersenne_twister[kk+1]&LOWER_MASK);
53 mersenne_twister[kk] = mersenne_twister[kk+(M-MERSENNE_TWISTER_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
54 }
55 y = (mersenne_twister[MERSENNE_TWISTER_N-1]&UPPER_MASK)|(mersenne_twister[0]&LOWER_MASK);
56 mersenne_twister[MERSENNE_TWISTER_N-1] = mersenne_twister[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
57
58 mersenne_twister_index = 0;
59 }
60
61
62 // returns current seed value
get_random_seed()63 uint32 get_random_seed()
64 {
65 if (mersenne_twister_index >= MERSENNE_TWISTER_N) { /* generate N words at one time */
66 MTgenerate();
67 }
68 return mersenne_twister[mersenne_twister_index];
69 }
70
71
72
73 /* generates a random number on [0,0xffffffff]-interval */
simrand_plain()74 uint32 simrand_plain()
75 {
76 uint32 y;
77
78 if (mersenne_twister_index >= MERSENNE_TWISTER_N) { /* generate N words at one time */
79 MTgenerate();
80 }
81 y = mersenne_twister[mersenne_twister_index++];
82
83 /* Tempering */
84 y ^= (y >> 11);
85 y ^= (y << 7) & 0x9d2c5680UL;
86 y ^= (y << 15) & 0xefc60000UL;
87 y ^= (y >> 18);
88
89 return y;
90 }
91
92
93 /* generates a random number on [0,max-1]-interval */
simrand(const uint32 max)94 uint32 simrand(const uint32 max)
95 {
96 assert( (random_origin&INTERACTIVE_RANDOM) == 0 );
97
98 if(max<=1) { // may rather assert this?
99 return 0;
100 }
101 return simrand_plain() % max;
102 }
103
104
clear_random_mode(uint16 mode)105 void clear_random_mode( uint16 mode )
106 {
107 random_origin &= ~mode;
108 }
109
110
set_random_mode(uint16 mode)111 void set_random_mode( uint16 mode )
112 {
113 random_origin |= mode;
114 }
115
116
get_random_mode()117 uint16 get_random_mode()
118 {
119 return random_origin;
120 }
121
122
123 static uint32 async_rand_seed = 12345678+dr_time();
124
125 // simpler simrand for anything not game critical (like UI)
sim_async_rand(uint32 max)126 uint32 sim_async_rand( uint32 max )
127 {
128 if( max==0 ) {
129 return 0;
130 }
131 async_rand_seed *= 3141592621u;
132 async_rand_seed ++;
133
134 return (async_rand_seed >> 8) % max;
135 }
136
137
138
139 static uint32 noise_seed = 0;
140
setsimrand(uint32 seed,uint32 ns)141 uint32 setsimrand(uint32 seed,uint32 ns)
142 {
143 uint32 old_noise_seed = noise_seed;
144
145 if(seed!=0xFFFFFFFF) {
146 init_genrand( seed );
147 async_rand_seed = seed+dr_time();
148 random_origin = 0;
149 }
150 if(noise_seed!=0xFFFFFFFF) {
151 noise_seed = ns*15731;
152 }
153 return old_noise_seed;
154 }
155
156
int_noise(const sint32 x,const sint32 y)157 static double int_noise(const sint32 x, const sint32 y)
158 {
159 uint32 n = (uint32)x + (uint32)y*101U + noise_seed;
160
161 n = (n<<13) ^ n;
162 return 1.0 - (double)((n * (n * n * 15731U + 789221U) + 1376312589U) & 0x7fffffff) / 1073741824.0;
163 }
164
165
166 static float *map = 0;
167 static sint32 map_w=0;
168
init_perlin_map(sint32 w,sint32 h)169 void init_perlin_map( sint32 w, sint32 h )
170 {
171 map_w = w+2;
172 map = new float[map_w*(h+2)];
173 for( sint32 y=0; y<h+2; y++ ) {
174 for( sint32 x=0; x<map_w; x++ ) {
175 map[x+(y*map_w)] = (float)int_noise( x-1, y-1 );
176 }
177 }
178 }
179
180
exit_perlin_map()181 void exit_perlin_map()
182 {
183 map_w = 0;
184 delete [] map;
185 map = 0;
186 }
187
188
189 #define map_noise(x,y) (0+map[(x)+1+((y)+1)*map_w])
190
191
smoothed_noise(const int x,const int y)192 static double smoothed_noise(const int x, const int y)
193 {
194 /* this gives a very smooth world */
195 if(map) {
196 const double corners =
197 map_noise(x-1, y-1)+map_noise(x+1, y-1)+map_noise(x-1, y+1)+map_noise(x+1, y+1);
198
199 const double sides =
200 map_noise(x-1, y) + map_noise(x+1, y) + map_noise(x, y-1) + map_noise(x, y+1);
201
202 const double center = map_noise(x, y);
203
204 return (corners + sides+sides + center*4.0) / 16.0;
205 }
206 else {
207 const double corners =
208 int_noise(x-1, y-1)+int_noise(x+1, y-1)+int_noise(x-1, y+1)+int_noise(x+1, y+1);
209
210 const double sides =
211 int_noise(x-1, y) + int_noise(x+1, y) + int_noise(x, y-1) + int_noise(x, y+1);
212
213 const double center = int_noise(x,y);
214
215 return (corners + sides+sides + center*4.0) / 16.0;
216 }
217
218
219 /* a hilly world
220 const double sides = ( int_noise(x-1, y) + int_noise(x+1, y) +
221 int_noise(x, y-1) + int_noise(x, y+1) );
222
223 const double center = int_noise(x, y);
224
225 return (sides+sides + center*4) / 8.0;
226 */
227
228 // this gives very hilly world
229 // return int_noise(x,y);
230 }
231
linear_interpolate(const double a,const double b,const double x)232 static double linear_interpolate(const double a, const double b, const double x)
233 {
234 // return a*(1.0-x) + b*x;
235 // return a - a*x + b*x;
236 return a + x*(b-a);
237 }
238
239
interpolated_noise(const double x,const double y)240 static double interpolated_noise(const double x, const double y)
241 {
242 // The function floor is needed because (int) rounds always towards zero,
243 // but we need integer_x be the biggest integer not bigger than x.
244 // So (int) -1.5 = -1.0
245 // But (int)floor(-1.5) = -2.0
246 // Modified 2008/10/17 by Gerd Wachsmuth
247 const int integer_X = (int)floor(x);
248 const int integer_Y = (int)floor(y);
249
250 const double fractional_X = x - (double)integer_X;
251 const double fractional_Y = y - (double)integer_Y;
252
253 const double v1 = smoothed_noise(integer_X, integer_Y);
254 const double v2 = smoothed_noise(integer_X + 1, integer_Y);
255 const double v3 = smoothed_noise(integer_X, integer_Y + 1);
256 const double v4 = smoothed_noise(integer_X + 1, integer_Y + 1);
257
258 const double i1 = linear_interpolate(v1 , v2 , fractional_X);
259 const double i2 = linear_interpolate(v3 , v4 , fractional_X);
260
261 return linear_interpolate(i1 , i2 , fractional_Y);
262 }
263
264
265 /**
266 * x,y Point coordinates
267 * p Persistence
268 */
perlin_noise_2D(const double x,const double y,const double p)269 double perlin_noise_2D(const double x, const double y, const double p)
270 {
271 double total = 0.0;
272 for( int i=0; i<6; i++ ) {
273 const double frequency = (double)(1 << i);
274 const double amplitude = pow(p, (double)i);
275 total += interpolated_noise( (x * frequency) / 64.0, (y * frequency) / 64.0) * amplitude;
276 }
277
278 return total;
279 }
280
281
282 // compute integer log10
log10(uint32 v)283 uint32 log10(uint32 v)
284 {
285 return ( (log2(v) + 1) * 1233) >> 12; // 1 / log_2(10) ~~ 1233 / 4096
286 }
287
288
log2(uint32 n)289 uint32 log2(uint32 n)
290 {
291 uint32 i = (n & 0xffff0000) ? 16 : 0;
292 if( (n >>= i) & 0xff00 ) {
293 i |= 8;
294 n >>= 8;
295 }
296 if( n & 0xf0 ) {
297 i |= 4;
298 n >>= 4;
299 }
300 if( n & 0xC ) {
301 i |= 2;
302 n >>= 2;
303 }
304 return i | (n >> 1);
305 }
306
307
308 // compute integer sqrt
sqrt_i32(uint32 num)309 uint32 sqrt_i32(uint32 num)
310 {
311 // taken from http://en.wikipedia.org/wiki/Methods_of_computing_square_roots
312 uint32 res = 0;
313 uint32 bit = 1 << 30; // The second-to-top bit is set: 1<<14 for short
314
315 // "bit" starts at the highest power of four <= the argument.
316 while( bit > num ) {
317 bit >>= 2;
318 }
319
320 while( bit != 0 ) {
321 if( num >= res + bit ) {
322 num -= res + bit;
323 res = (res >> 1) + bit;
324 }
325 else {
326 res >>= 1;
327 }
328 bit >>= 2;
329 }
330 return res;
331 }
332
333
334 // compute integer sqrt
sqrt_i64(uint64 num)335 uint64 sqrt_i64(uint64 num)
336 {
337 // taken from http://en.wikipedia.org/wiki/Methods_of_computing_square_roots
338 uint64 res = 0;
339 uint64 bit = (uint64)1 << 62; // The second-to-top bit is set: 1<<14 for short
340
341 // "bit" starts at the highest power of four <= the argument.
342 while( bit > num ) {
343 bit >>= 2;
344 }
345
346 while( bit != 0 ) {
347 if( num >= res + bit ) {
348 num -= res + bit;
349 res = (res >> 1) + bit;
350 }
351 else {
352 res >>= 1;
353 }
354 bit >>= 2;
355 }
356 return res;
357 }
358