small text and comment changes
[rocksndiamonds.git] / src / game_bd / bd_random.c
1 /* GLIB - Library of useful routines for C programming
2  * Copyright (C) 1995-1997  Peter Mattis, Spencer Kimball and Josh MacDonald
3  *
4  * SPDX-License-Identifier: LGPL-2.1-or-later
5  *
6  * This library is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Lesser General Public
8  * License as published by the Free Software Foundation; either
9  * version 2.1 of the License, or (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this library; if not, see <http://www.gnu.org/licenses/>.
18  */
19
20 /* Originally developed and coded by Makoto Matsumoto and Takuji
21  * Nishimura.  Please mail <matumoto@math.keio.ac.jp>, if you're using
22  * code from this file in your own programs or libraries.
23  * Further information on the Mersenne Twister can be found at
24  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
25  * This code was adapted to glib by Sebastian Wilhelmi.
26  */
27
28 /*
29  * Modified by the GLib Team and others 1997-2000.  See the AUTHORS
30  * file for a list of people on the GLib Team.  See the ChangeLog
31  * files for a list of changes.  These files are distributed with
32  * GLib at ftp://ftp.gtk.org/pub/gtk/.
33  */
34
35 #include <math.h>
36 #include <errno.h>
37 #include <stdio.h>
38 #include <string.h>
39 #include <sys/types.h>
40 #include <unistd.h>
41 #include <time.h>
42
43 #if defined(PLATFORM_WINDOWS)
44 #include <process.h>    // for getpid()
45 #endif
46
47 #include "main_bd.h"
48
49
50 /**
51  * GdRand:
52  *
53  * The GdRand struct is an opaque data structure. It should only be
54  * accessed through the gd_rand_* functions.
55  **/
56
57 // Period parameters
58 #define N 624
59 #define M 397
60 #define MATRIX_A   0x9908b0df // constant vector a
61 #define UPPER_MASK 0x80000000 // most significant w-r bits
62 #define LOWER_MASK 0x7fffffff // least significant r bits
63
64 // Tempering parameters
65 #define TEMPERING_MASK_B 0x9d2c5680
66 #define TEMPERING_MASK_C 0xefc60000
67 #define TEMPERING_SHIFT_U(y)  (y >> 11)
68 #define TEMPERING_SHIFT_S(y)  (y << 7)
69 #define TEMPERING_SHIFT_T(y)  (y << 15)
70 #define TEMPERING_SHIFT_L(y)  (y >> 18)
71
72 struct _GdRand
73 {
74   unsigned int mt[N]; // the array for the state vector
75   unsigned int mti;
76 };
77
78 /**
79  * gd_rand_new_with_seed: (constructor)
80  * @seed: a value to initialize the random number generator
81  *
82  * Creates a new random number generator initialized with @seed.
83  *
84  * Returns: (transfer full): the new #GdRand
85  **/
86 GdRand *
87 gd_rand_new_with_seed (unsigned int seed)
88 {
89   GdRand *rand = checked_calloc(sizeof(GdRand));
90   gd_rand_set_seed (rand, seed);
91   return rand;
92 }
93
94 /**
95  * gd_rand_new_with_seed_array: (constructor)
96  * @seed: an array of seeds to initialize the random number generator
97  * @seed_length: an array of seeds to initialize the random number
98  *     generator
99  *
100  * Creates a new random number generator initialized with @seed.
101  *
102  * Returns: (transfer full): the new #GdRand
103  *
104  * Since: 2.4
105  */
106 GdRand *
107 gd_rand_new_with_seed_array (const unsigned int *seed, unsigned int seed_length)
108 {
109   GdRand *rand = checked_calloc(sizeof(GdRand));
110   gd_rand_set_seed_array (rand, seed, seed_length);
111
112   return rand;
113 }
114
115 /**
116  * gd_rand_new: (constructor)
117  *
118  * Creates a new random number generator initialized with a seed taken
119  * either from `/dev/urandom` (if existing) or from the current time
120  * (as a fallback).
121  *
122  * On Windows, the seed is taken from rand_s().
123  *
124  * Returns: (transfer full): the new #GdRand
125  */
126 GdRand *
127 gd_rand_new (void)
128 {
129   unsigned int seed[4];
130
131 #if defined(PLATFORM_UNIX)
132   static boolean dev_urandom_exists = TRUE;
133
134   if (dev_urandom_exists)
135   {
136     FILE *dev_urandom;
137
138     do
139     {
140       dev_urandom = fopen ("/dev/urandom", "rbe");
141     }
142     while (dev_urandom == NULL && errno == EINTR);
143
144     if (dev_urandom)
145     {
146       int r;
147
148       setvbuf (dev_urandom, NULL, _IONBF, 0);
149
150       do
151       {
152         errno = 0;
153         r = fread (seed, sizeof (seed), 1, dev_urandom);
154       }
155       while (errno == EINTR);
156
157       if (r != 1)
158         dev_urandom_exists = FALSE;
159
160       fclose (dev_urandom);
161     }
162     else
163     {
164       dev_urandom_exists = FALSE;
165     }
166   }
167
168   if (!dev_urandom_exists)
169   {
170     struct timespec now;
171
172     clock_gettime(CLOCK_REALTIME, &now);
173
174     seed[0] = now.tv_sec;
175     seed[1] = now.tv_nsec;
176     seed[2] = getpid ();
177     seed[3] = getppid ();
178   }
179 #else // PLATFORM_WINDOWS
180   /* rand_s() is only available since Visual Studio 2005 and
181    * MinGW-w64 has a wrapper that will emulate rand_s() if it's not in msvcrt
182    */
183 #if (defined(_MSC_VER) && _MSC_VER >= 1400) || defined(__MINGW64_VERSION_MAJOR)
184   size_t i;
185
186   for (i = 0; i < ARRAY_SIZE(seed); i++)
187     rand_s (&seed[i]);
188 #else
189 #warning Using insecure seed for random number generation because of missing rand_s() in Windows XP
190   GTimeVal now;
191
192   gd_get_current_time (&now);
193
194   seed[0] = now.tv_sec;
195   seed[1] = now.tv_usec;
196   seed[2] = getpid ();
197   seed[3] = 0;
198 #endif
199
200 #endif
201
202   return gd_rand_new_with_seed_array (seed, 4);
203 }
204
205 /**
206  * gd_rand_free:
207  * @rand_: a #GdRand
208  *
209  * Frees the memory allocated for the #GdRand.
210  */
211 void
212 gd_rand_free (GdRand *rand)
213 {
214   if (rand != NULL)
215     checked_free(rand);
216 }
217
218 /**
219  * gd_rand_copy:
220  * @rand_: a #GdRand
221  *
222  * Copies a #GdRand into a new one with the same exact state as before.
223  * This way you can take a snapshot of the random number generator for
224  * replaying later.
225  *
226  * Returns: (transfer full): the new #GdRand
227  *
228  * Since: 2.4
229  */
230 GdRand *
231 gd_rand_copy (GdRand *rand)
232 {
233   GdRand *new_rand;
234
235   if (rand == NULL)
236     return NULL;
237
238   new_rand = checked_calloc(sizeof(GdRand));
239   memcpy(new_rand, rand, sizeof(GdRand));
240
241   return new_rand;
242 }
243
244 /**
245  * gd_rand_set_seed:
246  * @rand_: a #GdRand
247  * @seed: a value to reinitialize the random number generator
248  *
249  * Sets the seed for the random number generator #GdRand to @seed.
250  */
251 void
252 gd_rand_set_seed (GdRand *rand, unsigned int seed)
253 {
254   if (rand == NULL)
255     return;
256
257   // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
258   // In the previous version (see above), MSBs of the
259   // seed affect only MSBs of the array mt[].
260
261   rand->mt[0] = seed;
262   for (rand->mti = 1; rand->mti < N; rand->mti++)
263     rand->mt[rand->mti] = 1812433253UL *
264       (rand->mt[rand->mti - 1] ^ (rand->mt[rand->mti - 1] >> 30)) + rand->mti;
265 }
266
267 /**
268  * gd_rand_set_seed_array:
269  * @rand_: a #GdRand
270  * @seed: array to initialize with
271  * @seed_length: length of array
272  *
273  * Initializes the random number generator by an array of longs.
274  * Array can be of arbitrary size, though only the first 624 values
275  * are taken.  This function is useful if you have many low entropy
276  * seeds, or if you require more then 32 bits of actual entropy for
277  * your application.
278  *
279  * Since: 2.4
280  */
281 void
282 gd_rand_set_seed_array (GdRand *rand, const unsigned int *seed, unsigned int seed_length)
283 {
284   unsigned int i, j, k;
285
286   if (rand == NULL || seed_length < 1)
287     return;
288
289   gd_rand_set_seed (rand, 19650218UL);
290
291   i = 1;
292   j = 0;
293   k = (N > seed_length ? N : seed_length);
294
295   for (; k; k--)
296   {
297     rand->mt[i] = ((rand->mt[i] ^ ((rand->mt[i - 1] ^ (rand->mt[i - 1] >> 30)) * 1664525UL))
298                    + seed[j] + j);      // non linear
299     rand->mt[i] &= 0xffffffffUL;        // for WORDSIZE > 32 machines
300     i++;
301     j++;
302
303     if (i >= N)
304     {
305       rand->mt[0] = rand->mt[N - 1];
306       i = 1;
307     }
308
309     if (j >= seed_length)
310       j = 0;
311   }
312
313   for (k = N - 1; k; k--)
314   {
315     rand->mt[i] = ((rand->mt[i] ^ ((rand->mt[i - 1] ^ (rand->mt[i - 1] >> 30)) * 1566083941UL))
316                    - i);                // non linear
317     rand->mt[i] &= 0xffffffffUL;        // for WORDSIZE > 32 machines
318     i++;
319
320     if (i >= N)
321     {
322       rand->mt[0] = rand->mt[N - 1];
323       i = 1;
324     }
325   }
326
327   rand->mt[0] = 0x80000000UL;           // MSB is 1; assuring non-zero initial array
328 }
329
330 /**
331  * gd_rand_boolean:
332  * @rand_: a #GdRand
333  *
334  * Returns a random #boolean from @rand_.
335  * This corresponds to an unbiased coin toss.
336  *
337  * Returns: a random #boolean
338  */
339 /**
340  * gd_rand_int:
341  * @rand_: a #GdRand
342  *
343  * Returns the next random unsigned int from @rand_ equally distributed over
344  * the range [0..2^32-1].
345  *
346  * Returns: a random number
347  */
348 unsigned int
349 gd_rand_int (GdRand *rand)
350 {
351   unsigned int y;
352   static const unsigned int mag01[2] = { 0x0, MATRIX_A };
353   // mag01[x] = x * MATRIX_A  for x=0,1
354
355   if (rand == NULL)
356     return 0;
357
358   if (rand->mti >= N)
359   {
360     // generate N words at one time
361     int kk;
362
363     for (kk = 0; kk < N - M; kk++)
364     {
365       y = (rand->mt[kk] & UPPER_MASK) | (rand->mt[kk + 1] & LOWER_MASK);
366       rand->mt[kk] = rand->mt[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
367     }
368
369     for (; kk < N - 1; kk++)
370     {
371       y = (rand->mt[kk] & UPPER_MASK) | (rand->mt[kk + 1] & LOWER_MASK);
372       rand->mt[kk] = rand->mt[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
373     }
374
375     y = (rand->mt[N - 1] & UPPER_MASK) | (rand->mt[0] & LOWER_MASK);
376     rand->mt[N - 1] = rand->mt[M - 1] ^ (y >> 1) ^ mag01[y & 0x1];
377
378     rand->mti = 0;
379   }
380
381   y = rand->mt[rand->mti++];
382   y ^= TEMPERING_SHIFT_U(y);
383   y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
384   y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
385   y ^= TEMPERING_SHIFT_L(y);
386
387   return y;
388 }
389
390 // transform [0..2^32] -> [0..1]
391 #define GD_RAND_DOUBLE_TRANSFORM 2.3283064365386962890625e-10
392
393 /**
394  * gd_rand_int_range:
395  * @rand_: a #GdRand
396  * @begin: lower closed bound of the interval
397  * @end: upper open bound of the interval
398  *
399  * Returns the next random #int from @rand_ equally distributed over
400  * the range [@begin..@end-1].
401  *
402  * Returns: a random number
403  */
404 int
405 gd_rand_int_range (GdRand *rand, int begin, int end)
406 {
407   unsigned int dist = end - begin;
408   unsigned int random = 0;
409
410   if (rand == NULL || end <= begin)
411     return begin;
412
413   if (dist == 0)
414     return begin;
415
416   /* maxvalue is set to the predecessor of the greatest
417    * multiple of dist less or equal 2^32.
418    */
419   unsigned int maxvalue;
420   if (dist <= 0x80000000u) // 2^31
421   {
422     // maxvalue = 2^32 - 1 - (2^32 % dist)
423     unsigned int leftover = (0x80000000u % dist) * 2;
424     if (leftover >= dist) leftover -= dist;
425     maxvalue = 0xffffffffu - leftover;
426   }
427   else
428   {
429     maxvalue = dist - 1;
430   }
431
432   do
433     random = gd_rand_int (rand);
434   while (random > maxvalue);
435
436   random %= dist;
437
438   return begin + random;
439 }
440
441 /**
442  * gd_rand_double:
443  * @rand_: a #GRand
444  *
445  * Returns the next random #double from @rand_ equally distributed over
446  * the range [0..1).
447  *
448  * Returns: a random number
449  */
450 double
451 gd_rand_double (GdRand *rand)
452 {
453   /* We set all 52 bits after the point for this, not only the first
454      32. That's why we need two calls to gd_rand_int */
455   double retval = gd_rand_int(rand) * GD_RAND_DOUBLE_TRANSFORM;
456   retval = (retval + gd_rand_int(rand)) * GD_RAND_DOUBLE_TRANSFORM;
457
458   /* The following might happen due to very bad rounding luck, but
459    * actually this should be more than rare, we just try again then */
460   if (retval >= 1.0)
461     return gd_rand_double (rand);
462
463   return retval;
464 }
465
466 /**
467  * gd_rand_double_range:
468  * @rand_: a #GRand
469  * @begin: lower closed bound of the interval
470  * @end: upper open bound of the interval
471  *
472  * Returns the next random #double from @rand_ equally distributed over
473  * the range [@begin..@end).
474  *
475  * Returns: a random number
476  */
477 double
478 gd_rand_double_range (GdRand *rand, double begin, double end)
479 {
480   double r;
481
482   r = gd_rand_double(rand);
483
484   return r * end - (r - 1) * begin;
485 }
486
487 static GdRand *
488 get_global_random (void)
489 {
490   static GdRand *global_random;
491
492   // called while locked
493   if (!global_random)
494     global_random = gd_rand_new();
495
496   return global_random;
497 }
498
499 /**
500  * gd_random_boolean:
501  *
502  * Returns a random #gboolean.
503  * This corresponds to an unbiased coin toss.
504  *
505  * Returns: a random #gboolean
506  */
507 /**
508  * gd_random_int:
509  *
510  * Return a random #guint32 equally distributed over the range
511  * [0..2^32-1].
512  *
513  * Returns: a random number
514  */
515 unsigned int
516 gd_random_int (void)
517 {
518   unsigned int result;
519
520   result = gd_rand_int(get_global_random());
521
522   return result;
523 }
524
525 /**
526  * gd_random_int_range:
527  * @begin: lower closed bound of the interval
528  * @end: upper open bound of the interval
529  *
530  * Returns a random #int equally distributed over the range
531  * [@begin..@end-1].
532  *
533  * Returns: a random number
534  */
535 int
536 gd_random_int_range (int begin, int end)
537 {
538   int result;
539
540   result = gd_rand_int_range (get_global_random(), begin, end);
541
542   return result;
543 }
544
545 /**
546  * gd_random_double:
547  *
548  * Returns a random #double equally distributed over the range [0..1).
549  *
550  * Returns: a random number
551  */
552 double
553 gd_random_double (void)
554 {
555   double result;
556
557   result = gd_rand_double(get_global_random());
558
559   return result;
560 }
561
562 /**
563  * gd_random_double_range:
564  * @begin: lower closed bound of the interval
565  * @end: upper open bound of the interval
566  *
567  * Returns a random #double equally distributed over the range
568  * [@begin..@end).
569  *
570  * Returns: a random number
571  */
572 double
573 gd_random_double_range (double begin, double end)
574 {
575   double result;
576
577   result = gd_rand_double_range(get_global_random(), begin, end);
578
579   return result;
580 }