commit last state
[ppam-mpi.git] / code / src / Util / rng.c
1 // ISAAC PRNG
2 // Initial code by Bob Jenkins, March 1996 - modified in 2008
3 // Reference: http://burtleburtle.net/bob/rand/isaacafa.html
4 // Further slightly modified by Benjamin Auder, 2013
5
6 #include "Util/rng.h"
7 #include <stdlib.h>
8 #include <time.h>
9
10 // Internal state
11 static uint32_t mm[256];
12 static uint32_t aa, bb, cc;
13
14 // initialize the (pseudo-)random generator
15 void init_rng(int flag)
16 {
17 // 'bootstrap' isaac with basic bad PRNG
18 srand(time(NULL));
19
20 aa = rand() & 0xffffffff;
21 bb = rand() & 0xffffffff;
22 cc = rand() & 0xffffffff;
23
24 uint32_t a, b, c, d, e, f, g, h, i;
25 // The golden ratio:
26 a = b = c = d = e = f = g = h = 0x9e3779b9;
27 // Scramble it
28 for (i = 0; i < 4; ++i)
29 mix(a, b, c, d, e, f, g, h);
30
31 // Fill in mm[] with messy stuff (the seed)
32 uint32_t seedArray[256];
33 for (i = 0; i < 256; ++i) mm[i] = seedArray[i] = rand() & 0xffffffff;
34 for (i = 0; i < 256; i += 8)
35 {
36 if (flag)
37 {
38 // Use all the information in the seed
39 a += seedArray[i ];
40 b += seedArray[i + 1];
41 c += seedArray[i + 2];
42 d += seedArray[i + 3];
43 e += seedArray[i + 4];
44 f += seedArray[i + 5];
45 g += seedArray[i + 6];
46 h += seedArray[i + 7];
47 }
48 mix(a, b, c, d, e, f, g, h);
49 mm[i ] = a;
50 mm[i + 1] = b;
51 mm[i + 2] = c;
52 mm[i + 3] = d;
53 mm[i + 4] = e;
54 mm[i + 5] = f;
55 mm[i + 6] = g;
56 mm[i + 7] = h;
57 }
58
59 if (flag)
60 {
61 // Do a second pass to make all of the seed affect all of mm
62 for (i = 0; i < 256; i += 8)
63 {
64 a += mm[i ];
65 b += mm[i + 1];
66 c += mm[i + 2];
67 d += mm[i + 3];
68 e += mm[i + 4];
69 f += mm[i + 5];
70 g += mm[i + 6];
71 h += mm[i + 7];
72 mix(a, b, c, d, e, f, g, h);
73 mm[i ] = a;
74 mm[i + 1] = b;
75 mm[i + 2] = c;
76 mm[i + 3] = d;
77 mm[i + 4] = e;
78 mm[i + 5] = f;
79 mm[i + 6] = g;
80 mm[i + 7] = h;
81 }
82 }
83 }
84
85 // return a (pseudo-)random integer
86 uint32_t get_rand_int()
87 {
88 // TODO: register variables ? (x,y,i)
89 uint32_t x, y;
90 static uint32_t i = 0;
91
92 if (i == 0)
93 {
94 cc = cc + 1; // cc just gets incremented once per 256 results
95 bb = bb + cc; // then combined with bb
96 }
97
98 x = mm[i];
99 switch (i % 4)
100 {
101 case 0:
102 aa = aa^(aa << 13);
103 break;
104 case 1:
105 aa = aa^(aa >> 6);
106 break;
107 case 2:
108 aa = aa^(aa << 2);
109 break;
110 case 3:
111 aa = aa^(aa >> 16);
112 break;
113 }
114
115 // NOTE: bits 2..9 are chosen from x but 10..17 are chosen
116 // from y. The only important thing here is that 2..9 and 10..17
117 // don't overlap. 2..9 and 10..17 were then chosen for speed in
118 // the optimized version (rand.c) */
119 // See http://burtleburtle.net/bob/rand/isaac.html
120 // for further explanations and analysis.
121
122 aa = mm[(i + 128) % 256] + aa;
123 mm[i] = y = mm[(x >> 2) % 256] + aa + bb;
124 bb = mm[(y >> 10) % 256] + x;
125
126 i = (i + 1) % 256;
127
128 return bb;
129 }
130
131 // return a (pseudo-)random real number in [0,1]
132 Real get_rand_real()
133 {
134 return (Real) (INV_RANDMAX * (double) get_rand_int());
135 }