1 #ifdef __GNUG__
2 #pragma implementation
3 #endif
4 #include "ACG.h"
5 #include "assert.h"
6
7 //
8 // This is an extension of the older implementation of Algorithm M
9 // which I previously supplied. The main difference between this
10 // version and the old code are:
11 //
12 // + Andres searched high & low for good constants for
13 // the LCG.
14 //
15 // + theres more bit chopping going on.
16 //
17 // The following contains his comments.
18 //
19 // agn@UNH.CS.CMU.EDU sez..
20 //
21 // The generator below is based on 2 well known
22 // methods: Linear Congruential (LCGs) and Additive
23 // Congruential generators (ACGs).
24 //
25 // The LCG produces the longest possible sequence
26 // of 32 bit random numbers, each being unique in
27 // that sequence (it has only 32 bits of state).
28 // It suffers from 2 problems: a) Independence
29 // isnt great, that is the (n+1)th number is
30 // somewhat related to the preceding one, unlike
31 // flipping a coin where knowing the past outcomes
32 // dont help to predict the next result. b)
33 // Taking parts of a LCG generated number can be
34 // quite non-random: for example, looking at only
35 // the least significant byte gives a permuted
36 // 8-bit counter (that has a period length of only
37 // 256). The advantage of an LCA is that it is
38 // perfectly uniform when run for the entire period
39 // length (and very uniform for smaller sequences
40 // too, if the parameters are chosen carefully).
41 //
42 // ACGs have extremly long period lengths and
43 // provide good independence. Unfortunately,
44 // uniformity isnt not too great. Furthermore, I
45 // didnt find any theoretically analysis of ACGs
46 // that addresses uniformity.
47 //
48 // The RNG given below will return numbers
49 // generated by an LCA that are permuted under
50 // control of a ACG. 2 permutations take place: the
51 // 4 bytes of one LCG generated number are
52 // subjected to one of 16 permutations selected by
53 // 4 bits of the ACG. The permutation a such that
54 // byte of the result may come from each byte of
55 // the LCG number. This effectively destroys the
56 // structure within a word. Finally, the sequence
57 // of such numbers is permuted within a range of
58 // 256 numbers. This greatly improves independence.
59 //
60 //
61 // Algorithm M as describes in Knuths "Art of Computer Programming",
62 // Vol 2. 1969
63 // is used with a linear congruential generator (to get a good uniform
64 // distribution) that is permuted with a Fibonacci additive congruential
65 // generator to get good independence.
66 //
67 // Bit, byte, and word distributions were extensively tested and pass
68 // Chi-squared test near perfect scores (>7E8 numbers tested, Uniformity
69 // assumption holds with probability > 0.999)
70 //
71 // Run-up tests for on 7E8 numbers confirm independence with
72 // probability > 0.97.
73 //
74 // Plotting random points in 2d reveals no apparent structure.
75 //
76 // Autocorrelation on sequences of 5E5 numbers (A(i) = SUM X(n)*X(n-i),
77 // i=1..512)
78 // results in no obvious structure (A(i) ~ const).
79 //
80 // Except for speed and memory requirements, this generator outperforms
81 // random() for all tests. (random() scored rather low on uniformity tests,
82 // while independence test differences were less dramatic).
83 //
84 // AGN would like to..
85 // thanks to M.Mauldin, H.Walker, J.Saxe and M.Molloy for inspiration & help.
86 //
87 // And I would (DGC) would like to thank Donald Kunth for AGN for letting me
88 // use his extensions in this implementation.
89 //
90
91 //
92 // Part of the table on page 28 of Knuth, vol II. This allows us
93 // to adjust the size of the table at the expense of shorter sequences.
94 //
95
96 static randomStateTable[][3] = {
97 {3,7,16}, {4,9, 32}, {3,10, 32}, {1,11, 32}, {1,15,64}, {3,17,128},
98 {7,18,128}, {3,20,128}, {2,21, 128}, {1,22, 128}, {5,23, 128}, {3,25, 128},
99 {2,29, 128}, {3,31, 128}, {13,33, 256}, {2,35, 256}, {11,36, 256},
100 {14,39,256}, {3,41,256}, {9,49,256}, {3,52,256}, {24,55,256}, {7,57, 256},
101 {19,58,256}, {38,89,512}, {17,95,512}, {6,97,512}, {11,98,512}, {-1,-1,-1} };
102
103 //
104 // spatial permutation table
105 // RANDOM_PERM_SIZE must be a power of two
106 //
107
108 #define RANDOM_PERM_SIZE 64
109 unsigned long randomPermutations[RANDOM_PERM_SIZE] = {
110 0xffffffff, 0x00000000, 0x00000000, 0x00000000, // 3210
111 0x0000ffff, 0x00ff0000, 0x00000000, 0xff000000, // 2310
112 0xff0000ff, 0x0000ff00, 0x00000000, 0x00ff0000, // 3120
113 0x00ff00ff, 0x00000000, 0xff00ff00, 0x00000000, // 1230
114
115 0xffff0000, 0x000000ff, 0x00000000, 0x0000ff00, // 3201
116 0x00000000, 0x00ff00ff, 0x00000000, 0xff00ff00, // 2301
117 0xff000000, 0x00000000, 0x000000ff, 0x00ffff00, // 3102
118 0x00000000, 0x00000000, 0x00000000, 0xffffffff, // 2103
119
120 0xff00ff00, 0x00000000, 0x00ff00ff, 0x00000000, // 3012
121 0x0000ff00, 0x00000000, 0x00ff0000, 0xff0000ff, // 2013
122 0x00000000, 0x00000000, 0xffffffff, 0x00000000, // 1032
123 0x00000000, 0x0000ff00, 0xffff0000, 0x000000ff, // 1023
124
125 0x00000000, 0xffffffff, 0x00000000, 0x00000000, // 0321
126 0x00ffff00, 0xff000000, 0x00000000, 0x000000ff, // 0213
127 0x00000000, 0xff000000, 0x0000ffff, 0x00ff0000, // 0132
128 0x00000000, 0xff00ff00, 0x00000000, 0x00ff00ff // 0123
129 };
130
131 //
132 // SEED_TABLE_SIZE must be a power of 2
133 //
134 #define SEED_TABLE_SIZE 32
135 static unsigned long seedTable[SEED_TABLE_SIZE] = {
136 0xbdcc47e5, 0x54aea45d, 0xec0df859, 0xda84637b,
137 0xc8c6cb4f, 0x35574b01, 0x28260b7d, 0x0d07fdbf,
138 0x9faaeeb0, 0x613dd169, 0x5ce2d818, 0x85b9e706,
139 0xab2469db, 0xda02b0dc, 0x45c60d6e, 0xffe49d10,
140 0x7224fea3, 0xf9684fc9, 0xfc7ee074, 0x326ce92a,
141 0x366d13b5, 0x17aaa731, 0xeb83a675, 0x7781cb32,
142 0x4ec7c92d, 0x7f187521, 0x2cf346b4, 0xad13310f,
143 0xb89cff2b, 0x12164de1, 0xa865168d, 0x32b56cdf
144 };
145
146 //
147 // The LCG used to scramble the ACG
148 //
149 //
150 // LC-parameter selection follows recommendations in
151 // "Handbook of Mathematical Functions" by Abramowitz & Stegun 10th, edi.
152 //
153 // LC_A = 251^2, ~= sqrt(2^32) = 66049
154 // LC_C = result of a long trial & error series = 3907864577
155 //
156
157 static const unsigned long LC_A = 66049;
158 static const unsigned long LC_C = 3907864577;
LCG(unsigned long x)159 static inline unsigned long LCG(unsigned long x)
160 {
161 return( x * LC_A + LC_C );
162 }
163
164
ACG(unsigned long seed,int size)165 ACG::ACG(unsigned long seed, int size)
166 {
167
168 initialSeed = seed;
169
170 //
171 // Determine the size of the state table
172 //
173
174 for (register int l = 0;
175 randomStateTable[l][0] != -1 && randomStateTable[l][1] < size;
176 l++);
177
178 if (randomStateTable[l][1] == -1) {
179 l--;
180 }
181
182 initialTableEntry = l;
183
184 stateSize = randomStateTable[ initialTableEntry ][ 1 ];
185 auxSize = randomStateTable[ initialTableEntry ][ 2 ];
186
187 //
188 // Allocate the state table & the auxillary table in a single malloc
189 //
190
191 state = new unsigned long[stateSize + auxSize];
192 auxState = &state[stateSize];
193
194 reset();
195 }
196
197 //
198 // Initialize the state
199 //
200 void
reset()201 ACG::reset()
202 {
203 register unsigned long u;
204
205 if (initialSeed > -1 && initialSeed < SEED_TABLE_SIZE) {
206 u = seedTable[ initialSeed ];
207 } else {
208 u = initialSeed ^ seedTable[ initialSeed & (SEED_TABLE_SIZE-1) ];
209 }
210
211
212 j = randomStateTable[ initialTableEntry ][ 0 ] - 1;
213 k = randomStateTable[ initialTableEntry ][ 1 ] - 1;
214
215 register int i;
216 for(i = 0; i < stateSize; i++) {
217 state[i] = u = LCG(u);
218 }
219
220 for (i = 0; i < auxSize; i++) {
221 auxState[i] = u = LCG(u);
222 }
223
224 k = u % stateSize;
225 int tailBehind = (stateSize - randomStateTable[ initialTableEntry ][ 0 ]);
226 j = k - tailBehind;
227 if (j < 0) {
228 j += stateSize;
229 }
230
231 lcgRecurr = u;
232
233 assert(sizeof(double) == 2 * sizeof(long));
234 }
235
~ACG()236 ACG::~ACG()
237 {
238 if (state) delete state;
239 state = 0;
240 // don't delete auxState, it's really an alias for state.
241 }
242
243 //
244 // Returns 32 bits of random information.
245 //
246
asLong()247 unsigned long ACG::asLong()
248 {
249 unsigned long result = state[k] + state[j];
250 state[k] = result;
251 j = (j <= 0) ? (stateSize-1) : (j-1);
252 k = (k <= 0) ? (stateSize-1) : (k-1);
253
254 short int auxIndex = (result >> 24) & (auxSize - 1);
255 register unsigned long auxACG = auxState[auxIndex];
256 auxState[auxIndex] = lcgRecurr = LCG(lcgRecurr);
257
258 //
259 // 3c is a magic number. We are doing four masks here, so we
260 // do not want to run off the end of the permutation table.
261 // This insures that we have always got four entries left.
262 //
263 register unsigned long *perm = & randomPermutations[result & 0x3c];
264
265 result = *(perm++) & auxACG;
266 result |= *(perm++) & ((auxACG << 24)
267 | ((auxACG >> 8)& 0xffffff));
268 result |= *(perm++) & ((auxACG << 16)
269 | ((auxACG >> 16) & 0xffff));
270 result |= *(perm++) & ((auxACG << 8)
271 | ((auxACG >> 24) & 0xff));
272
273 return(result);
274 }
275