1 /* rng.c
2 This file contains the code for a high-quality random number
3 generator written by Don Knuth. The auxilliary routine
4 ran_arr_cycle() has been modified slightly, and ran_init() is new.
5
6 To use it:
7
8 0. #include "rng.h" (or "naututil.h" if you are using nauty)
9
10 1. Call ran_init(seed), where seed is any long integer.
11 This step is optional, but if you don't use it you
12 will always get the same sequence of random numbers.
13
14 2. For each random number, use the NEXTRAN macro. It will
15 give a random value in the range 0..2^30-1. Alternatively,
16 KRAN(k) will have a random value in the range 0..k-1.
17 KRAN(k) actually gives you NEXTRAN mod k, so it is not
18 totally uniform if k is very large. In that case, you
19 can use the slightly slower GETKRAN(k,var) to set the
20 variable var to a better random number from 0..k-1.
21
22 Brendan McKay, July 2002. Fixed Nov 2002 on advice of DEK.
23
24 */
25
26 /* This program by D E Knuth is in the public domain and freely copyable
27 * AS LONG AS YOU MAKE ABSOLUTELY NO CHANGES!
28 * It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
29 * (or in the errata to the 2nd edition --- see
30 * http://www-cs-faculty.stanford.edu/~knuth/taocp.html
31 * in the changes to Volume 2 on pages 171 and following). */
32
33 /* N.B. The MODIFICATIONS introduced in the 9th printing (2002) are
34 included here; there's no backwards compatibility with the original. */
35
36 /* If you find any bugs, please report them immediately to
37 * taocp@cs.stanford.edu
38 * (and you will be rewarded if the bug is genuine). Thanks! */
39
40 /************ see the book for explanations and caveats! *******************/
41 /************ in particular, you need two's complement arithmetic **********/
42
43 #define KK 100 /* the long lag */
44 #define LL 37 /* the short lag */
45 #define MM (1L<<30) /* the modulus */
46 #define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
47
48 long ran_x[KK]; /* the generator state */
49
50 #ifdef __STDC__
ran_array(long aa[],int n)51 void ran_array(long aa[],int n)
52 #else
53 void ran_array(aa,n) /* put n new random numbers in aa */
54 long *aa; /* destination */
55 int n; /* array length (must be at least KK) */
56 #endif
57 {
58 int i,j;
59 for (j=0;j<KK;j++) aa[j]=ran_x[j];
60 for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
61 for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
62 for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
63 }
64
65 /* the following routines are from exercise 3.6--15 */
66 /* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */
67
68 #define QUALITY 1009 /* recommended quality level for high-res use */
69 static long ran_arr_buf[QUALITY];
70 static long ran_arr_dummy=-1, ran_arr_started=-1;
71 long *ran_arr_ptr=&ran_arr_dummy; /* the next random number, or -1 */
72
73 #define TT 70 /* guaranteed separation between streams */
74 #define is_odd(x) ((x)&1) /* units bit of x */
75
76 #ifdef __STDC__
ran_start(long seed)77 void ran_start(long seed)
78 #else
79 void ran_start(seed) /* do this before using ran_array */
80 long seed; /* selector for different streams */
81 #endif
82 {
83 int t,j;
84 long x[KK+KK-1]; /* the preparation buffer */
85 long ss=(seed+2)&(MM-2);
86 for (j=0;j<KK;j++) {
87 x[j]=ss; /* bootstrap the buffer */
88 ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
89 }
90 x[1]++; /* make x[1] (and only x[1]) odd */
91 for (ss=seed&(MM-1),t=TT-1; t; ) {
92 for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */
93 for (j=KK+KK-2;j>=KK;j--)
94 x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]),
95 x[j-KK]=mod_diff(x[j-KK],x[j]);
96 if (is_odd(ss)) { /* "multiply by z" */
97 for (j=KK;j>0;j--) x[j]=x[j-1];
98 x[0]=x[KK]; /* shift the buffer cyclically */
99 x[LL]=mod_diff(x[LL],x[KK]);
100 }
101 if (ss) ss>>=1; else t--;
102 }
103 for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
104 for (;j<KK;j++) ran_x[j-LL]=x[j];
105 for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
106 ran_arr_ptr=&ran_arr_started;
107 }
108
109 void
ran_init(long seed)110 ran_init(long seed) /* Added by BDM: use instead of ran_start. */
111 /* But this is less important with this version */
112 {
113 ran_start((unsigned long)seed % (MM-2));
114 }
115
116 #define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
117
118 long
ran_arr_cycle(void)119 ran_arr_cycle(void)
120 /* Modified by BDM to automatically initialise
121 if no explicit initialisation has been done */
122 {
123 if (ran_arr_ptr==&ran_arr_dummy)
124 ran_start(314159L); /* the user forgot to initialize */
125
126 ran_array(ran_arr_buf,QUALITY);
127
128 ran_arr_buf[KK]=-1;
129 ran_arr_ptr=ran_arr_buf+1;
130 return ran_arr_buf[0];
131 }
132