1 /*
2 * rng_uvag.c
3 *
4 * Copyright(c) 2004 by Alex Hay - zenjew@hotmail.com
5 *
6 * From:
7 *
8 * Private communication, 6/5/07.
9 *
10 * GSL packaged for dieharder by Robert G. Brown 6/5/07
11 *
12 * Alex has agree to make this available under the GSL so it can be
13 * integrated with the dieharder package.
14 *
15 * Here is pretty much the only documentation on this generator in
16 * existence, at least until dieharder-based studies enable a real
17 * publication to occur. The following is from Alex:
18 *
19 * UVAG: The Universal Virtual Array Generator
20 *
21 * This generator is universal in the sense that a single algorithm
22 * can be used to produce random numbers with any integer data TYPE.
23 * Thus it is trivial to extend it to 64 bits and beyond. If and when
24 * long long = 128 bits, UVAG is ready to go.
25 *
26 * The array is virtual in that when TYPE>char the 256 inner TYPES overlap
27 * producing an effectively huge "virtual" inner complexity in a very small
28 * physical memory. When TYPE = uint, every cycle changes 4-7 (depending on
29 * where rp is pointing in the array) virtual TYPEs in the array. When
30 * TYPE = long long each cycle changes 8-15 virtual TYPES at a physical
31 * memory cost of just 4 bytes over uint for the entire array.
32 *
33 * As sindex is incremented, its next potential value has been changed many
34 * times since it was last encountered along with the entire inner state.
35 * Not very rigorous but then I'm no mathematician. UVAG is something I
36 * would dub a chaotic prng. It doesn't yield easily to mathematical
37 * analysis unlike well researched prngs like MT and RC4. I have
38 * experimented using fewer than 256 seeds. Fewer than 100 do well in
39 * DIEHARD. The cost is in speed due to an extra mod.
40 *
41 * Where did UVAG come from?
42 *
43 * It evolved over years of playing around. In its previous incarnation, rp
44 * was moved randomly and proved (as Makoto Matsumoto predicted) to zero out
45 * after extended cycling. Systematically stepping thru the array with sindex
46 * seems to eliminate the problem although I have not consulted further with
47 * Dr Matsumoto.
48 *
49 *========================================================================
50 *
51 * This program is free software; you can redistribute it and/or modify
52 * it under the terms of the GNU General Public License as published by
53 * the Free Software Foundation; either version 2 of the License, or (at
54 * your option) any later version.
55 *
56 * This program is distributed in the hope that it will be useful, but
57 * WITHOUT ANY WARRANTY; without even the implied warranty of
58 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
59 * General Public License for more details.
60 *
61 * You should have received a copy of the GNU General Public License
62 * along with this program; if not, write to the Free Software
63 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
64 */
65
66 #include <dieharder/libdieharder.h>
67
68 /*
69 * This is a wrapping of the uvag rng
70 */
71
72 static unsigned long int uvag_get (void *vstate);
73 static double uvag_get_double (void *vstate);
74 static void uvag_set (void *vstate, unsigned long int s);
75
76 typedef struct
77 {
78 int i;
79 }
80 uvag_state_t;
81
82 /*
83 * UVAG specific defines.
84 *
85 * UVAG provides globally as *rndint, an effectively infinite sequence of
86 * TYPEs, uniformly distributed (0, 2**(8*sizeof(TYPE))-1). We need to
87 * return uints in our GSL-compatible wrapper, at least at the moment.
88 */
89
90 /* #define TYPE long long || char || short || int || long */
91 #define TYPE uint
92 #define WORD sizeof(TYPE)
93
94 /*
95 * Global variables and data for UVAG
96 */
97
98 TYPE *rp;
99 TYPE rndint;
100 unsigned char sindex, svec[255 + WORD]; /* 256 overlapping TYPE seeds */
101
102 static inline unsigned long int
uvag_get(void * vstate)103 uvag_get (void *vstate)
104 {
105
106 /*
107 * Returns a 32-bit unsigned integer produced by the UVAG
108 */
109 rp = (TYPE *)(svec + svec[sindex+=1]);
110 rndint += *rp += rndint;
111 return rndint;
112
113 }
114
115 static double
uvag_get_double(void * vstate)116 uvag_get_double (void *vstate)
117 {
118 return uvag_get (vstate) / (double) UINT_MAX;
119 }
120
uvag_set(void * vstate,unsigned long int s)121 static void uvag_set (void *vstate, unsigned long int s) {
122
123 /* Initialize automaton using specified seed. */
124 uvag_state_t *state __attribute__((unused)) = (uvag_state_t *) vstate;
125
126 uint i, array_len = 255 + WORD, tot, seed_seed, tmp8;
127 unsigned char key[256], *kp, temp;
128 gsl_rng *seed_rng; /* random number generator used to seed uvag */
129
130 /*
131 * Preload the array with 1-byte integers
132 */
133 for (i=0; i<array_len; i++){
134 svec[i]=i;
135 }
136
137 /*
138 * OK, here we have to modify Alex's algorithm. The GSL requires a
139 * single seed, unsigned long int in type. Alex requires a key string
140 * 256 characters long (that is, 64 uints long). We therefore have to
141 * bootstrap off of an existing, deterministic GSL RNG to fill the seed
142 * string from a single permitted seed. Note that type 12 is the
143 * mt19937_1999 generator, basically one of the best in the world -- not
144 * that it matters.
145 */
146 seed_rng = gsl_rng_alloc(dh_rng_types[14]);
147 seed_seed = s;
148 gsl_rng_set(seed_rng,seed_seed);
149 random_max = gsl_rng_max(seed_rng);
150 rmax = random_max;
151 rmax_bits = 0;
152 rmax_mask = 0;
153 while(rmax){
154 rmax >>= 1;
155 rmax_mask = rmax_mask << 1;
156 rmax_mask++;
157 rmax_bits++;
158 }
159 for(i=0;i<256;i++){
160 /* if(i%32 == 0) printf("\n"); */
161 get_rand_bits(&tmp8,sizeof(uint),8,seed_rng);
162 if(i!=255){
163 key[i] = tmp8;
164 } else {
165 key[i] = 0;
166 }
167 /* printf("%02x",key[i]); */
168 }
169 /* printf("\n"); */
170
171 kp = key;
172 tot = 0;
173 for(i=0; i<array_len; i++) { /* shuffle seeds */
174 tot += *kp++;
175 tot%=array_len;
176 temp = svec[tot];
177 svec[tot] = svec[i];
178 svec[i] = temp;
179 /* wrap around key[] */
180 if(*kp == '\0') {
181 kp = key;
182 }
183 }
184 /* For debugging of the original load'n'shuffle
185 printf("svec = ");
186 for(i=0;i<256;i++){
187 if(i%32 == 0) printf("\n");
188 printf("%02x|",svec[i]);
189 }
190 printf("\n");
191 */
192
193 sindex = 0;
194 rndint = 0;
195
196 }
197
198 static const gsl_rng_type uvag_type =
199 {"uvag", /* name */
200 UINT_MAX, /* RAND_MAX */
201 0, /* RAND_MIN */
202 sizeof (uvag_state_t),
203 &uvag_set,
204 &uvag_get,
205 &uvag_get_double};
206
207 const gsl_rng_type *gsl_rng_uvag = &uvag_type;
208