// ISAAC PRNG // Initial code by Bob Jenkins, March 1996 - modified in 2008 // Reference: http://burtleburtle.net/bob/rand/isaacafa.html // Further slightly modified by Benjamin Auder, 2013 #include "Util/rng.h" #include #include // Internal state static uint32_t mm[256]; static uint32_t aa, bb, cc; // initialize the (pseudo-)random generator void init_rng(int flag) { // 'bootstrap' isaac with basic bad PRNG srand(time(NULL)); aa = rand() & 0xffffffff; bb = rand() & 0xffffffff; cc = rand() & 0xffffffff; uint32_t a, b, c, d, e, f, g, h, i; // The golden ratio: a = b = c = d = e = f = g = h = 0x9e3779b9; // Scramble it for (i = 0; i < 4; ++i) mix(a, b, c, d, e, f, g, h); // Fill in mm[] with messy stuff (the seed) uint32_t seedArray[256]; for (i = 0; i < 256; ++i) mm[i] = seedArray[i] = rand() & 0xffffffff; for (i = 0; i < 256; i += 8) { if (flag) { // Use all the information in the seed a += seedArray[i ]; b += seedArray[i + 1]; c += seedArray[i + 2]; d += seedArray[i + 3]; e += seedArray[i + 4]; f += seedArray[i + 5]; g += seedArray[i + 6]; h += seedArray[i + 7]; } mix(a, b, c, d, e, f, g, h); mm[i ] = a; mm[i + 1] = b; mm[i + 2] = c; mm[i + 3] = d; mm[i + 4] = e; mm[i + 5] = f; mm[i + 6] = g; mm[i + 7] = h; } if (flag) { // Do a second pass to make all of the seed affect all of mm for (i = 0; i < 256; i += 8) { a += mm[i ]; b += mm[i + 1]; c += mm[i + 2]; d += mm[i + 3]; e += mm[i + 4]; f += mm[i + 5]; g += mm[i + 6]; h += mm[i + 7]; mix(a, b, c, d, e, f, g, h); mm[i ] = a; mm[i + 1] = b; mm[i + 2] = c; mm[i + 3] = d; mm[i + 4] = e; mm[i + 5] = f; mm[i + 6] = g; mm[i + 7] = h; } } } // return a (pseudo-)random integer uint32_t get_rand_int() { // TODO: register variables ? (x,y,i) uint32_t x, y; static uint32_t i = 0; if (i == 0) { cc = cc + 1; // cc just gets incremented once per 256 results bb = bb + cc; // then combined with bb } x = mm[i]; switch (i % 4) { case 0: aa = aa^(aa << 13); break; case 1: aa = aa^(aa >> 6); break; case 2: aa = aa^(aa << 2); break; case 3: aa = aa^(aa >> 16); break; } // NOTE: bits 2..9 are chosen from x but 10..17 are chosen // from y. The only important thing here is that 2..9 and 10..17 // don't overlap. 2..9 and 10..17 were then chosen for speed in // the optimized version (rand.c) */ // See http://burtleburtle.net/bob/rand/isaac.html // for further explanations and analysis. aa = mm[(i + 128) % 256] + aa; mm[i] = y = mm[(x >> 2) % 256] + aa + bb; bb = mm[(y >> 10) % 256] + x; i = (i + 1) % 256; return bb; } // return a (pseudo-)random real number in [0,1] Real get_rand_real() { return (Real) (INV_RANDMAX * (double) get_rand_int()); }