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