1 /**********************************************************************
2 random_util.c
3 **********************************************************************
4
5 random_util - Pseudo-random number generation routines.
6 Copyright ©2000-2004, Stewart Adcock <stewart@linux-domain.com>
7 All rights reserved.
8
9 The latest version of this program should be available at:
10 http://gaul.sourceforge.net/
11
12 This program is free software; you can redistribute it and/or modify
13 it under the terms of the GNU General Public License as published by
14 the Free Software Foundation; either version 2 of the License, or
15 (at your option) any later version. Alternatively, if your project
16 is incompatible with the GPL, I will probably agree to requests
17 for permission to use the terms of any other license.
18
19 This program is distributed in the hope that it will be useful, but
20 WITHOUT ANY WARRANTY WHATSOEVER.
21
22 A full copy of the GNU General Public License should be in the file
23 "COPYING" provided with this distribution; if not, see:
24 http://www.gnu.org/
25
26 **********************************************************************
27
28 Synopsis: Portable routines for generating pseudo random numbers.
29
30 Why use these instead of the system routines?
31 (a) Useful interface functions. (can specify ranges or
32 specific distributions)
33 (b) System independant. (Generated numbers are
34 reproducible across system types.)
35 (c) Enables saving and restoring state.
36
37 SLang intrinsic function wrappers are provided if the
38 HAVE_SLANG constant is set to 1.
39
40 The algorithm I selected is the Mitchell and Moore
41 variant of the standard additive number generator.
42 The required array of numbers is populated by a
43 using single seed value by using a linear
44 congruential pseudo random number generator.
45
46 This routines have been tested and provide the same
47 output (within the limits of computational precision)
48 on (at least) the following platforms:
49
50 o Intel x86 (Linux, OpenBSD, FreeBSD)
51 o AMD x86-64 (Linux, FreeBSD)
52 o IBM Power3 (AIX)
53 o IBM PowerPC (Linux)
54 o MIPS R4K, R10K, R12K (IRIX)
55 o Alpha EV56 (Tru64), EV7 (Linux)
56 o SPARC Ultra-4 (Solaris)
57
58 This code should be thread safe.
59
60 For OpenMP code, USE_OPENMP must be defined and
61 random_init() or random_seed() must be called BY A
62 SINGLE THREAD prior to any other function.
63
64 References: Standard additive number generator and the linear
65 congruential algorithm:
66 Knuth D.E., "The Art of Computer Programming",
67 Vol 2, 2nd ed. (1998)
68
69 General information:
70 Press W., Flannery B.P., Teukolsky S.A., Vetterling W.T.,
71 "Numerical Recipes in C: The Art of Scientific Computing"
72 Cambridge University Press, 2nd ed. (1992)
73
74 Usage: o First call random_init().
75
76 o Call random_seed() to seed the PRNG (like srand()).
77 Alternatively, random_tseed() to use system clock for seed.
78
79 o random_rand() is a rand() replacement, which is available
80 for use but I would recommend the wrapper functions:
81 random_int(), random_int_range() etc.
82
83 o random_get_state() and random_set_state() may be used
84 to set, save, restore, and query the current state.
85
86 These functions can be tested by compiling with
87 something like:
88 gcc -o testrand random_util.c -DRANDOM_UTIL_TEST
89
90 To do: Gaussian/Boltzmann distributions.
91 Need a proper test of randomness.
92 Properly implement stack of states.
93 I could do with some handy documentation.
94
95 **********************************************************************/
96
97 #include "gaul/random_util.h"
98
99 /*
100 * PRNG constants. Don't mess with these values willie-nillie.
101 */
102 #define RANDOM_MM_ALPHA 55
103 #define RANDOM_MM_BETA 24
104 #define RANDOM_LC_ALPHA 3
105 #define RANDOM_LC_BETA 257
106 #define RANDOM_LC_GAMMA RANDOM_RAND_MAX
107
108 /*
109 * Global state variable stack.
110 * (Implemented using singly-linked list.)
111 */
112 /*
113 static GSList *state_list=NULL;
114 */
115 static random_state current_state;
116 static boolean is_initialised=FALSE;
117
118 THREAD_LOCK_DEFINE_STATIC(random_state_lock);
119
120 /**********************************************************************
121 random_rand()
122 Synopsis: Replacement for the standard rand().
123 Returns a new pseudo-random value from the sequence, in
124 the range 0 to RANDOM_RAND_MAX inclusive, and updates
125 global state for next call. size should be non-zero,
126 and state should be initialized.
127 parameters:
128 return: none
129 last updated: 30/12/00
130 **********************************************************************/
131
random_rand(void)132 unsigned int random_rand(void)
133 {
134 unsigned int val;
135
136 if (!is_initialised) die("Neither random_init() or random_seed() have been called.");
137
138 THREAD_LOCK(random_state_lock);
139
140 val = (current_state.v[current_state.j]+current_state.v[current_state.k])
141 & RANDOM_RAND_MAX;
142
143 current_state.x = (current_state.x+1) % RANDOM_NUM_STATE_VALS;
144 current_state.j = (current_state.j+1) % RANDOM_NUM_STATE_VALS;
145 current_state.k = (current_state.k+1) % RANDOM_NUM_STATE_VALS;
146 current_state.v[current_state.x] = val;
147
148 THREAD_UNLOCK(random_state_lock);
149
150 return val;
151 }
152
153
154 /**********************************************************************
155 random_seed()
156 synopsis: Set seed for pseudo random number generator.
157 Uses the linear congruential algorithm to fill
158 state array.
159 parameters: const unsigned int seed Seed value.
160 return: none
161 last updated: 04 May 2004
162 **********************************************************************/
163
random_seed(const unsigned int seed)164 void random_seed(const unsigned int seed)
165 {
166 int i;
167
168 #if USE_OPENMP == 1
169 if (is_initialised == FALSE)
170 {
171 omp_init_lock(&random_state_lock);
172 is_initialised = TRUE;
173 }
174 #else
175 is_initialised = TRUE;
176 #endif
177
178 THREAD_LOCK(random_state_lock);
179
180 current_state.v[0]=(seed & RANDOM_RAND_MAX);
181
182 for(i=1; i<RANDOM_NUM_STATE_VALS; i++)
183 current_state.v[i] = (RANDOM_LC_ALPHA * current_state.v[i-1]
184 + RANDOM_LC_BETA) & RANDOM_RAND_MAX;
185
186 current_state.j = 0;
187 current_state.k = RANDOM_MM_ALPHA-RANDOM_MM_BETA;
188 current_state.x = RANDOM_MM_ALPHA-0;
189
190 THREAD_UNLOCK(random_state_lock);
191
192 return;
193 }
194
195
196 /**********************************************************************
197 random_tseed()
198 synopsis: Set seed for pseudo random number generator from
199 the system time.
200 parameters:
201 return: none
202 last updated: 28/01/01
203 **********************************************************************/
204
random_tseed(void)205 void random_tseed(void)
206 {
207 random_seed((unsigned int) (time(NULL)%RANDOM_RAND_MAX));
208
209 return;
210 }
211
212
213 /*************************************************************************
214 random_get_state()
215 synopsis: Retrieve current state. This can be used for saving
216 the current state.
217 parameters: none
218 return: random_state state
219 last updated: 16/05/00
220 *************************************************************************/
221
random_get_state(void)222 random_state random_get_state(void)
223 {
224 return current_state;
225 }
226
227
228 /*************************************************************************
229 random_set_state()
230 synopsis: Replace current state with specified state.
231 This can be used for restoring a saved state.
232 parameters: random_state state
233 return: none
234 last updated: 16/05/00
235 *************************************************************************/
236
random_set_state(random_state state)237 void random_set_state(random_state state)
238 {
239 current_state = state;
240
241 return;
242 }
243
244
245 /*************************************************************************
246 random_get_state_str()
247 synopsis: Retrieve current state encoded as a static string.
248 This can be used for saving the current state.
249 parameters: none
250 return: char *
251 last updated: 28/01/01
252 *************************************************************************/
253
random_get_state_str(void)254 char *random_get_state_str(void)
255 {
256 return (char *) ¤t_state;
257 }
258
259
260 /*************************************************************************
261 random_get_state_str_len()
262 synopsis: Retrieve the length of the string encoded current state.
263 This can be used for saving the current state.
264 parameters: none
265 return: char *
266 last updated: 28/01/01
267 *************************************************************************/
268
random_get_state_str_len(void)269 unsigned int random_get_state_str_len(void)
270 {
271 return (unsigned int) sizeof(random_state);
272 }
273
274
275 /*************************************************************************
276 random_set_state_str()
277 synopsis: Replace current state with specified state.
278 This can be used for restoring a saved state.
279 parameters: char *
280 return: none
281 last updated: 28/01/01
282 *************************************************************************/
283
random_set_state_str(char * state)284 void random_set_state_str(char *state)
285 {
286 /* This causes potential unaligned trap on Alpha CPUs. */
287 current_state = *((random_state *)state);
288
289 return;
290 }
291
292
293 /**********************************************************************
294 random_init()
295 synopsis: Initialise random number generators. Should be
296 called prior to any of these functions being used.
297 Seeds the PRNG with a seed of 1.
298 parameters: none
299 return: none
300 last updated: 30 May 2002
301 **********************************************************************/
302
random_init(void)303 void random_init(void)
304 {
305 random_seed(1);
306
307 #if RANDOM_DEBUG>0
308 printf("DEBUG: Random number routines initialised.\n");
309 #endif
310
311 return;
312 }
313
314
315 /**********************************************************************
316 random_isinit()
317 synopsis: Whether these routines have been initialised.
318 parameters:
319 return: none
320 last updated: 30/12/00
321 **********************************************************************/
322
random_isinit(void)323 boolean random_isinit(void)
324 {
325 return is_initialised;
326 }
327
328
329 /**********************************************************************
330 PRNG interface routines.
331 **********************************************************************/
332
333 /**********************************************************************
334 random_boolean()
335 synopsis: Return TRUE 50% of the time.
336 parameters:
337 return: none
338 last updated: 16/05/00
339 **********************************************************************/
340
random_boolean(void)341 boolean random_boolean(void)
342 {
343 return (random_rand() <= RANDOM_RAND_MAX/2);
344 }
345
346
347 /**********************************************************************
348 random_boolean_prob()
349 synopsis: Return TRUE (prob*100)% of the time.
350 parameters:
351 return: TRUE or FALSE.
352 last updated: 16/05/00
353 **********************************************************************/
354
random_boolean_prob(const double prob)355 boolean random_boolean_prob(const double prob)
356 {
357 return (random_rand() <= (unsigned int)(prob*(double)RANDOM_RAND_MAX));
358 }
359
360
361 /**********************************************************************
362 random_int()
363 synopsis: Return a random integer between 0 and (N-1) inclusive.
364 parameters:
365 return:
366 last updated: 05 Jun 2003
367 **********************************************************************/
368
random_int(const unsigned int max)369 unsigned int random_int(const unsigned int max)
370 {
371 /*
372 return (int)((double)random_rand()*max/RANDOM_RAND_MAX);
373 */
374 if (max==0)
375 return 0;
376 else
377 return random_rand()%max;
378 }
379
380
381 /**********************************************************************
382 random_int_range()
383 synopsis: Return a random integer between min and max-1 inclusive.
384 parameters:
385 return:
386 last updated: 05 Jun 2003
387 **********************************************************************/
388
random_int_range(const int min,const int max)389 int random_int_range(const int min, const int max)
390 {
391 /*
392 return (random_rand()*(max-min)/RANDOM_RAND_MAX + min);
393 */
394 if (max-min==0)
395 return max;
396 else
397 return min + (random_rand()%(max-min));
398 }
399
400
401 /**********************************************************************
402 random_double_full()
403 synopsis: Return a random double within the allowed range.
404 parameters:
405 return:
406 last updated: 16/06/01
407 **********************************************************************/
408
random_double_full(void)409 double random_double_full(void)
410 {
411 return ( ((double)random_rand()/(double)RANDOM_RAND_MAX)*
412 (DBL_MAX-DBL_MIN)+DBL_MIN );
413 }
414
415
416 /**********************************************************************
417 random_double()
418 synopsis: Return a random double within the specified range.
419 parameters:
420 return:
421 last updated: 22/01/01
422 **********************************************************************/
423
random_double(const double max)424 double random_double(const double max)
425 {
426 return ( max*(((double)random_rand())/(double)RANDOM_RAND_MAX) );
427 }
428
429
430 /**********************************************************************
431 random_double_range()
432 synopsis: Return a random double within the specified range.
433 parameters:
434 return:
435 last updated: 07/06/00
436 **********************************************************************/
437
random_double_range(const double min,const double max)438 double random_double_range(const double min, const double max)
439 {
440 return ( (max-min)*(((double)random_rand())/(double)RANDOM_RAND_MAX) + min );
441 }
442
443
444 /**********************************************************************
445 random_double_1()
446 synopsis: Return a random double within the range -1.0=>r>1.0
447 parameters:
448 return:
449 last updated: 21/02/01
450 **********************************************************************/
451
random_double_1(void)452 double random_double_1(void)
453 {
454 return ( 2.0*(((double)random_rand())/(double)RANDOM_RAND_MAX) - 1.0 );
455 }
456
457
458 /**********************************************************************
459 random_float_full()
460 synopsis: Return a random float within the allowed range.
461 parameters:
462 return:
463 last updated: 4 Dec 2001
464 **********************************************************************/
465
random_float_full(void)466 float random_float_full(void)
467 {
468 return ( ((float)random_rand()/(float)RANDOM_RAND_MAX)*
469 (DBL_MAX-DBL_MIN)+DBL_MIN );
470 }
471
472
473 /**********************************************************************
474 random_float()
475 synopsis: Return a random float within the specified range.
476 parameters:
477 return:
478 last updated: 4 Dec 2001
479 **********************************************************************/
480
random_float(const float max)481 float random_float(const float max)
482 {
483 return ( max*(((float)random_rand())/(float)RANDOM_RAND_MAX) );
484 }
485
486
487 /**********************************************************************
488 random_float_range()
489 synopsis: Return a random float within the specified range.
490 parameters:
491 return:
492 last updated: 4 Dec 2001
493 **********************************************************************/
494
random_float_range(const float min,const float max)495 float random_float_range(const float min, const float max)
496 {
497 return ( (max-min)*(((float)random_rand())/(float)RANDOM_RAND_MAX) + min );
498 }
499
500
501 /**********************************************************************
502 random_float_1()
503 synopsis: Return a random float within the range -1.0=>r>1.0
504 parameters:
505 return:
506 last updated: 4 Dec 2001
507 **********************************************************************/
508
random_float_1(void)509 float random_float_1(void)
510 {
511 return ( 2.0*(((float)random_rand())/(float)RANDOM_RAND_MAX) - 1.0 );
512 }
513
514
515 /**********************************************************************
516 random_float_unit_uniform()
517 synopsis: Return a pseudo-random number with a uniform
518 distribution in the range 0.0=>r>1.0
519 parameters:
520 return:
521 last updated: 4 Dec 2001
522 **********************************************************************/
523
random_float_unit_uniform(void)524 float random_float_unit_uniform(void)
525 {
526 return ( (((float)random_rand())/(float)RANDOM_RAND_MAX) );
527 }
528
529
530 /**********************************************************************
531 random_unit_uniform()
532 synopsis: Return a pseudo-random number with a uniform
533 distribution in the range 0.0=>r>1.0
534 parameters:
535 return:
536 last updated: 11/01/01
537 **********************************************************************/
538
random_unit_uniform(void)539 double random_unit_uniform(void)
540 {
541 return ( (((double)random_rand())/(double)RANDOM_RAND_MAX) );
542 }
543
544
545 /**********************************************************************
546 random_float_gaussian()
547 synopsis: Return a pseudo-random number with a normal
548 distribution with a given mean and standard devaiation.
549 parameters:
550 return:
551 last updated: 4 Dec 2001
552 **********************************************************************/
553
random_float_gaussian(const float mean,const float stddev)554 float random_float_gaussian(const float mean, const float stddev)
555 {
556 float q,u,v,x,y;
557
558 /*
559 * Generate P = (u,v) uniform in rectangular acceptance region.
560 */
561 do
562 {
563 u = 1.0-random_float_unit_uniform(); /* in range 0.0>u>=1.0 */
564 v = 1.7156 * (0.5 - random_float_unit_uniform()); /* in range 0.8578>v>=0.8578 */
565
566 /* Evaluate the quadratic form. */
567 x = u - 0.449871;
568 y = fabs(v) + 0.386595;
569 q = x * x + y * (0.19600 * y - 0.25472 * x);
570
571 /*
572 * Accept P if inside inner ellipse.
573 * Reject P if outside outer ellipse, or outside acceptance region.
574 */
575 } while ((q >= 0.27597) && ((q > 0.27846) || (v * v > -4.0 * log(u) * u * u)));
576
577 /* Return ratio of P's coordinates as the normal deviate. */
578 return (mean + 2.0 * stddev * v / u); /* I'm not entirely sure why this *2.0 factor is needed! */
579 }
580
581
582 /**********************************************************************
583 random_float_unit_gaussian()
584 synopsis: Random number with normal distribution, average 0.0,
585 deviation 1.0
586 parameters:
587 return:
588 last updated: 07 Jan 2002
589 **********************************************************************/
590
random_float_unit_gaussian(void)591 float random_float_unit_gaussian(void)
592 {
593 float r, u, v, fac;
594 static boolean set = FALSE;
595 static float dset;
596
597 if (set)
598 {
599 set = FALSE;
600 return dset;
601 }
602
603 do
604 {
605 u = 2.0 * random_float_unit_uniform() - 1.0;
606 v = 2.0 * random_float_unit_uniform() - 1.0;
607 r = u*u + v*v;
608 } while (r >= 1.0);
609
610 fac = sqrt(-2.0 * log(r) / r);
611 dset = v*fac;
612
613 return u*fac;
614 }
615
616 /* The original (thread-safe) version was this: */
617 #if 0
618 float random_float_unit_gaussian(void)
619 {
620 float r, u, v;
621
622 do
623 {
624 u = 2.0 * random_float_unit_uniform() - 1.0;
625 v = 2.0 * random_float_unit_uniform() - 1.0;
626 r = u*u + v*v;
627 } while (r >= 1.0);
628
629 return u*sqrt(-2.0 * log(r) / r);
630 }
631 #endif
632
633
634 /**********************************************************************
635 random_float_cauchy()
636 synopsis: Random number with a Cauchy/Lorentzian distribution.
637 parameters:
638 return:
639 last updated: 4 Dec 2001
640 **********************************************************************/
641
random_float_cauchy(void)642 float random_float_cauchy(void)
643 {
644 return tan(random_float_range(-PI/2,PI/2));
645 }
646
647
648 /**********************************************************************
649 random_float_exponential()
650 synopsis: Random number with an exponential distribution, mean
651 of 1.0.
652 parameters:
653 return:
654 last updated: 4 Dec 2001
655 **********************************************************************/
656
random_float_exponential(void)657 float random_float_exponential(void)
658 {
659 return -log(random_float_unit_uniform());
660 }
661
662
663 /**********************************************************************
664 random_gaussian()
665 synopsis: Return a pseudo-random number with a normal
666 distribution with a given mean and standard devaiation.
667 parameters:
668 return:
669 last updated: 11/01/01
670 **********************************************************************/
671
672 /*
673 Kinda based on: (But optimised quite a bit)
674
675 ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
676 THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
677 VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
678 The algorithm uses the ratio of uniforms method of A.J. Kinderman
679 and J.F. Monahan augmented with quadratic bounding curves.
680 */
random_gaussian(const double mean,const double stddev)681 double random_gaussian(const double mean, const double stddev)
682 {
683 double q,u,v,x,y;
684
685 /*
686 * Generate P = (u,v) uniform in rectangular acceptance region.
687 */
688 do
689 {
690 u = 1.0-random_unit_uniform(); /* in range 0.0>u>=1.0 */
691 v = 1.7156 * (0.5 - random_unit_uniform()); /* in range 0.8578>v>=0.8578 */
692
693 /* Evaluate the quadratic form. */
694 x = u - 0.449871;
695 y = fabs(v) + 0.386595;
696 q = x * x + y * (0.19600 * y - 0.25472 * x);
697
698 /*
699 * Accept P if inside inner ellipse.
700 * Reject P if outside outer ellipse, or outside acceptance region.
701 */
702 } while ((q >= 0.27597) && ((q > 0.27846) || (v * v > -4.0 * log(u) * u * u)));
703
704 /* Return ratio of P's coordinates as the normal deviate. */
705 return (mean + 2.0 * stddev * v / u); /* I'm not entirely sure why this *2.0 factor is needed! */
706 }
707
708
709 #if 0
710 /* Random number with normal distribution, average 0, deviation 1.
711 From Numerical Recipes. */
712 double random_unit_gaussian(void)
713 {
714 static boolean set = FALSE;
715 static double dset;
716 double fac, r, u, v;
717
718 if (!set)
719 {
720 do
721 {
722 u = 2.0 * random_unit_uniform() - 1.0;
723 v = 2.0 * random_unit_uniform() - 1.0;
724 r = u*u + v*v;
725 } while (r >= 1.0);
726
727 fac = sqrt(-2.0 * log(r) / r);
728 dset = u * fac;
729 set = TRUE;
730
731 return v * fac;
732 }
733 else
734 {
735 set = FALSE;
736 return dset;
737 }
738 }
739 #endif
740
741
742 /**********************************************************************
743 random_unit_gaussian()
744 synopsis: Random number with normal distribution, average 0.0,
745 deviation 1.0
746 parameters:
747 return:
748 last updated: 07 Jan 2002
749 **********************************************************************/
750
random_unit_gaussian(void)751 double random_unit_gaussian(void)
752 {
753 double r, u, v, fac;
754 static boolean set = FALSE;
755 static double dset;
756
757 if (set)
758 {
759 set = FALSE;
760 return dset;
761 }
762
763 do
764 {
765 u = 2.0 * random_unit_uniform() - 1.0;
766 v = 2.0 * random_unit_uniform() - 1.0;
767 r = u*u + v*v;
768 } while (r >= 1.0);
769
770 fac = sqrt(-2.0 * log(r) / r);
771 dset = v*fac;
772
773 return u*fac;
774 }
775
776 /* The original (thread-safe) version was this: */
777 #if 0
778 double random_unit_gaussian(void)
779 {
780 double r, u, v;
781
782 do
783 {
784 u = 2.0 * random_unit_uniform() - 1.0;
785 v = 2.0 * random_unit_uniform() - 1.0;
786 r = u*u + v*v;
787 } while (r >= 1.0);
788
789 return u*sqrt(-2.0 * log(r) / r);
790 }
791 #endif
792
793
794 /**********************************************************************
795 random_cauchy()
796 synopsis: Random number with a Cauchy/Lorentzian distribution.
797 parameters: none
798 return: double Random value.
799 last updated: 30/04/01
800 **********************************************************************/
801
random_cauchy(void)802 double random_cauchy(void)
803 {
804 return tan(random_double_range(-PI/2,PI/2));
805 }
806
807
808 /**********************************************************************
809 random_exponential()
810 synopsis: Random number with an exponential distribution, mean
811 of 1.0.
812 parameters: none
813 return: double Random value.
814 last updated: 30/04/01
815 **********************************************************************/
816
random_exponential(void)817 double random_exponential(void)
818 {
819 return -log(random_unit_uniform());
820 }
821
822
823 /**********************************************************************
824 random_int_permutation()
825 synopsis: Randomize an array of integers.
826 parameters: int size
827 int *iarray Source array.
828 int *oarray Destination array.
829 return: none
830 last updated: 14 Mar 2002
831 **********************************************************************/
832
random_int_permutation(const int size,int * iarray,int * oarray)833 void random_int_permutation(const int size, int *iarray, int *oarray)
834 {
835 int i,j=0; /* Loop variables over arrays. */
836 int pos; /* Randomly selected index. */
837
838 if (!iarray || !oarray) die("NULL pointer to int array passed.");
839
840 for (i=size-1; i>0; i--)
841 {
842 pos = random_int(i);
843 oarray[j++] = iarray[pos];
844 iarray[pos] = iarray[i];
845 }
846
847 oarray[j] = iarray[0];
848
849 return;
850 }
851
852
853 /**********************************************************************
854 random_diagnostics()
855 synopsis: Diagnostics.
856 parameters: none
857 return: none
858 last updated: 13 Mar 2002
859 **********************************************************************/
860
random_diagnostics(void)861 void random_diagnostics(void)
862 {
863 int i; /* Loop over PRNG array. */
864
865 printf("=== PRNG routines diagnostic information =====================\n");
866 printf("Version: %s\n", GA_VERSION_STRING);
867 printf("Build date: %s\n", GA_BUILD_DATE_STRING);
868 printf("Compilation machine characteristics:\n%s\n", GA_UNAME_STRING);
869 printf("--------------------------------------------------------------\n");
870 printf("RANDOM_DEBUG: %d\n", RANDOM_DEBUG);
871 printf("RANDOM_RAND_MAX: %u\n", RANDOM_RAND_MAX);
872 printf("RANDOM_NUM_STATE_VALS: %d\n", RANDOM_NUM_STATE_VALS);
873 #if HAVE_SLANG==1
874 printf("HAVE_SLANG: TRUE\n");
875 #else
876 printf("HAVE_SLANG: FALSE\n");
877 #endif
878 printf("--------------------------------------------------------------\n");
879 printf("structure sizeof\n");
880 printf("random_state: %lu\n", (unsigned long) sizeof(random_state));
881 printf("--------------------------------------------------------------\n");
882
883 if (is_initialised)
884 {
885 printf("Current state\n");
886 printf("j: %d k: %d x: %d v[%d]:\n",
887 current_state.j, current_state.k, current_state.x, RANDOM_NUM_STATE_VALS);
888 for (i=0; i<RANDOM_NUM_STATE_VALS; i++) printf("%u ", current_state.v[i]);
889 printf("\n");
890 }
891 else
892 {
893 printf("Not initialised.\n");
894 }
895
896 printf("==============================================================\n");
897
898 return;
899 }
900
901
902 /**********************************************************************
903 random_test()
904 synopsis: Testing.
905 parameters:
906 return: TRUE if all tests successful.
907 last updated: 30/12/00
908 **********************************************************************/
909
910 #define NUM_BINS 200
911 #define NUM_SAMPLES 1000000
912 #define NUM_CHISQ 20
913
random_test(void)914 boolean random_test(void)
915 {
916 unsigned int i, j, k; /* Loop variables. */
917 double r; /* Pseudo-random number. */
918 long bins[NUM_BINS]; /* Bin. */
919 double sum=0,sumsq=0; /* Stats. */
920 int numtrue=0, numfalse=0; /* Count booleans. */
921 unsigned int rchi=100; /* Number of bins in chisq test. */
922 unsigned int nchi=1000; /* Number of samples in chisq test. */
923 double chisq; /* Chisq error. */
924 double elimit = 2*sqrt((double)rchi); /* Chisq error limit. */
925 double nchi_rchi = (double)nchi / (double)rchi; /* Speed calculation. */
926 FILE *rfile=NULL; /* Handle for file of random integers. */
927
928 random_init();
929
930 printf("Testing random numbers.\n");
931
932 /*
933 * Uniform Distribution.
934 */
935 printf("Uniform distribution. Mean should be about 0.5.\n");
936
937 for (i=0;i<NUM_BINS;i++) bins[i] = 0;
938
939 for (i=0;i<NUM_SAMPLES;i++)
940 {
941 r = random_unit_uniform();
942 if (r >= 0.0 && r < 1.0)
943 {
944 bins[(int)(r*NUM_BINS)]++;
945 sum += r;
946 sumsq += SQU(r);
947 }
948 else
949 {
950 printf("Number generated out of range 0.0<=r<1.0.\n");
951 }
952 }
953 printf("Mean = %f\n", sum / NUM_SAMPLES);
954 printf("Standard deviation = %f\n", (sumsq - sum*sum/NUM_SAMPLES)/NUM_SAMPLES);
955
956 for (i=0;i<NUM_BINS;i++)
957 printf("%5.3f %ld\n", i/(double)NUM_BINS, bins[i]);
958
959 /*
960 * Gaussian Distribution.
961 */
962 printf("Gaussian distribution. Mean should be about 0.45. Standard deviation should be about 0.05.\n");
963
964 sum=0;
965 sumsq=0;
966
967 for (i=0;i<NUM_BINS;i++) bins[i] = 0;
968
969 for (i=0;i<NUM_SAMPLES;i++)
970 {
971 r = random_gaussian(0.45,0.05);
972 if (r >= 0.0 && r < 1.0)
973 {
974 bins[(int)(r*NUM_BINS)]++;
975 sum += r;
976 sumsq += SQU(r);
977 }
978 else
979 {
980 printf("Number generated out of range 0.0<=r<1.0.\n");
981 }
982 }
983 printf("Mean = %f\n", sum / NUM_SAMPLES);
984 printf("Standard deviation = %f\n", (sumsq - sum*sum/NUM_SAMPLES)/NUM_SAMPLES);
985
986 for (i=0;i<NUM_BINS;i++)
987 printf("%5.3f %ld\n", i/(double)NUM_BINS, bins[i]);
988
989 /*
990 * Unit Gaussian Distribution.
991 */
992 printf("Gaussian distribution. Mean should be about 0.0. Standard deviation should be about 1.0.\n");
993
994 sum=0;
995 sumsq=0;
996
997 for (i=0;i<NUM_BINS;i++) bins[i] = 0;
998
999 for (i=0;i<NUM_SAMPLES;i++)
1000 {
1001 r = random_unit_gaussian();
1002 if (r >= -5.0 && r < 5.0)
1003 {
1004 bins[(int)((r+5.0)*NUM_BINS)/10]++;
1005 sum += r;
1006 sumsq += SQU(r);
1007 }
1008 else
1009 {
1010 printf("Number generated out of range -5.0<=r<5.0.\n");
1011 }
1012 }
1013 printf("Mean = %f\n", sum / NUM_SAMPLES);
1014 printf("Standard deviation = %f\n", (sumsq - sum*sum/NUM_SAMPLES)/NUM_SAMPLES);
1015
1016 for (i=0;i<NUM_BINS;i++)
1017 printf("%5.3f %ld\n", -5.0+10*i/(double)NUM_BINS, bins[i]);
1018
1019 /*
1020 * Random Boolean.
1021 */
1022 printf("Random Booleans. Two counts should be approximately equal.\n");
1023
1024 for (i=0;i<NUM_SAMPLES;i++)
1025 {
1026 if ( random_boolean() )
1027 numtrue++;
1028 else
1029 numfalse++;
1030 }
1031 printf("TRUE/FALSE = %d/%d\n", numtrue, numfalse);
1032
1033 /*
1034 * Random int.
1035 */
1036 printf("Random Integers. The distribution should be approximately uniform.\n");
1037
1038 for (i=0;i<NUM_BINS;i++) bins[i] = 0;
1039
1040 for (i=0;i<NUM_SAMPLES;i++)
1041 bins[random_int(NUM_BINS)]++;
1042
1043 for (i=0;i<NUM_BINS;i++)
1044 printf("%u %ld\n", i, bins[i]);
1045
1046 /*
1047 * Chi squared test. This is the standard basic test for randomness of a PRNG.
1048 * We would expect any moethod to fail about about one out of ten runs.
1049 * The error is r*t/N - N and should be within 2*sqrt(r) of r.
1050 */
1051
1052 printf("Chi Squared Test of Random Integers. We would expect a couple of failures.\n");
1053
1054 if (rchi>NUM_BINS) die("Internal error: too few bins.");
1055
1056 for (j=0;j<NUM_CHISQ;j++)
1057 {
1058 printf("Run %u. chisq should be within %f of %u.\n", j, elimit, rchi);
1059 for(k=0; k<10; k++)
1060 {
1061 memset(bins, 0, rchi*sizeof(long));
1062 chisq = 0.0;
1063
1064 for(i=0; i<nchi; i++)
1065 bins[random_int(rchi)]++;
1066
1067 for(i=0; i<rchi; i++)
1068 chisq += SQU((double)bins[i] - nchi_rchi);
1069
1070 chisq /= nchi_rchi;
1071
1072 printf("chisq = %f - %s\n", chisq, fabs(chisq - rchi)>elimit?"FAILED":"PASSED");
1073 }
1074 }
1075
1076 printf("Creating file (\"randtest.dat\") of 5000 random integer numbers for external analysis.\n");
1077
1078 rfile = fopen("randtest.dat", "w");
1079
1080 for (i=0; i<5000; i++)
1081 {
1082 r = (double) random_rand();
1083 fprintf(rfile, "%f %f\n",
1084 /*i, r,*/
1085 (double)i/(double)5000, r/(double)RANDOM_RAND_MAX);
1086 }
1087
1088 fclose(rfile);
1089
1090 return TRUE;
1091 }
1092
1093
1094 /**********************************************************************
1095 SLang intrinsic wrapper functions.
1096
1097 FIXME: Problem with unsigned vs. signed ints.
1098 **********************************************************************/
1099 #if HAVE_SLANG==1
1100
1101 /* These function don't need wrappers:
1102 void random_tseed(void)
1103 void random_init(void)
1104 boolean random_is_init(void)
1105 boolean random_boolean(void)
1106 double random_unit_uniform(void)
1107 double random_unit_gaussian(void)
1108 void random_diagnostics(void)
1109 boolean random_test(void)
1110 */
1111
random_rand_wrapper(void)1112 int random_rand_wrapper(void)
1113 {
1114 return (int) random_rand();
1115 }
1116
1117
random_seed_wrapper(int * seed)1118 void random_seed_wrapper(int *seed)
1119 {
1120 random_seed((unsigned int) *seed);
1121 return;
1122 }
1123
1124
1125 #if 0
1126 THREAD_LOCK_DEFINE_STATIC(state_table_lock);
1127 TableStruct *state_table=NULL; /* Table for state handles. */
1128
1129 int random_get_state_wrapper(void)
1130 {
1131 int stateid; /* State handle. */
1132
1133 THREAD_LOCK(state_table_lock);
1134 if (!state_table) state_table=table_new();
1135
1136 stateid = table_add(state_table, (vpointer) current_state); /* FIXME: Need to copy state. */
1137 THREAD_UNLOCK(state_table_lock);
1138
1139 return stateid;
1140 }
1141
1142
1143 void random_set_state_wrapper(int stateid)
1144 {
1145 random_state state; /* State to restore. */
1146
1147 THREAD_LOCK(state_table_lock);
1148 state = table_get_data(state_table, stateid);
1149 THREAD_UNLOCK(state_table_lock);
1150
1151 if (!state) die("Unmatched state handle.");
1152
1153 random_set_state(state);
1154
1155 return;
1156 }
1157 #endif
1158
1159
random_boolean_prob_wrapper(double * prob)1160 boolean random_boolean_prob_wrapper(double *prob)
1161 {
1162 return random_boolean_prob(*prob);
1163 }
1164
1165
random_int_wrapper(int * max)1166 int random_int_wrapper(int *max)
1167 {
1168 return (int) random_int((unsigned int) *max);
1169 }
1170
1171
random_int_range_wrapper(int * min,int * max)1172 int random_int_range_wrapper(int *min, int *max)
1173 {
1174 return random_int_range(*min, *max);
1175 }
1176
1177
random_double_wrapper(double * max)1178 double random_double_wrapper(double *max)
1179 {
1180 return random_double(*max);
1181 }
1182
1183
random_double_range_wrapper(double * min,double * max)1184 double random_double_range_wrapper(double *min, double *max)
1185 {
1186 return random_double_range(*min, *max);
1187 }
1188
1189
random_gaussian_wrapper(double * mean,double * stddev)1190 double random_gaussian_wrapper(double *mean, double *stddev)
1191 {
1192 return random_gaussian(*mean, *stddev);
1193 }
1194
1195 #endif
1196