isaac.h revision 12855:588919e0e4aa 
1#ifndef __ISAAC_HPP
2#define __ISAAC_HPP
3
4
5/*
6
7    C++ TEMPLATE VERSION OF Robert J. Jenkins Jr.'s
8    ISAAC Random Number Generator.
9
10    Ported from vanilla C to to template C++ class
11    by Quinn Tyler Jackson on 16-23 July 1998.
12
13        quinn@qtj.net
14
15    The function for the expected period of this
16    random number generator, according to Jenkins is:
17
18        f(a,b) = 2**((a+b*(3+2^^a)-1)
19
20        (where a is ALPHA and b is bitwidth)
21
22    So, for a bitwidth of 32 and an ALPHA of 8,
23    the expected period of ISAAC is:
24
25        2^^(8+32*(3+2^^8)-1) = 2^^8295
26
27    Jackson has been able to run implementations
28    with an ALPHA as high as 16, or
29
30        2^^2097263
31
32*/
33
34
35typedef unsigned int UINT32;
36const UINT32 GOLDEN_RATIO = UINT32(0x9e3779b9);
37
38
39template <UINT32 ALPHA = (8)>
40class QTIsaac
41{
42  public:
43
44    typedef unsigned char byte;
45
46    struct randctx
47    {
48        randctx(void)
49        {
50            randrsl = new UINT32[N];
51            randmem = new UINT32[N];
52        }
53
54        ~randctx(void)
55        {
56            delete [] randrsl;
57            delete [] randmem;
58        }
59
60        UINT32 randcnt;
61        UINT32* randrsl;
62        UINT32* randmem;
63        UINT32 randa;
64        UINT32 randb;
65        UINT32 randc;
66	};
67
68    QTIsaac(UINT32 a = 0, UINT32 b = 0, UINT32 c = 0);
69    virtual ~QTIsaac(void);
70
71    UINT32 rand(void);
72    virtual void randinit(randctx* ctx, bool bUseSeed);
73    virtual void srand(
74		UINT32 a = 0, UINT32 b = 0, UINT32 c = 0, UINT32* s = NULL);
75
76    enum {N = (1<<ALPHA)};
77
78  protected:
79
80     virtual void isaac(randctx* ctx);
81
82     UINT32 ind(UINT32* mm, UINT32 x);
83     void rngstep(
84		UINT32 mix, UINT32& a, UINT32& b, UINT32*& mm, UINT32*& m,
85		UINT32*& m2, UINT32*& r, UINT32& x, UINT32& y);
86     virtual void shuffle(
87		UINT32& a, UINT32& b, UINT32& c, UINT32& d, UINT32& e, UINT32& f,
88		UINT32& g, UINT32& h);
89
90  private:
91    randctx m_rc;
92};
93
94
95template<UINT32 ALPHA>
96QTIsaac<ALPHA>::QTIsaac(UINT32 a, UINT32 b, UINT32 c) : m_rc()
97{
98    srand(a, b, c);
99}
100
101
102template<UINT32 ALPHA>
103QTIsaac<ALPHA>::~QTIsaac(void)
104{
105    // DO NOTHING
106}
107
108
109template<UINT32 ALPHA>
110void QTIsaac<ALPHA>::srand(UINT32 a, UINT32 b, UINT32 c, UINT32* s)
111{
112	for(int i = 0; i < N; i++)
113	{
114		m_rc.randrsl[i] = s != NULL ? s[i] : 0;
115	}
116
117	m_rc.randa = a;
118	m_rc.randb = b;
119	m_rc.randc = c;
120
121	randinit(&m_rc, true);
122}
123
124
125template<UINT32 ALPHA>
126inline UINT32 QTIsaac<ALPHA>::rand(void)
127{
128	return 0x7fffffff & (!m_rc.randcnt-- ?
129		(isaac(&m_rc), m_rc.randcnt=(N-1), m_rc.randrsl[m_rc.randcnt]) :
130		m_rc.randrsl[m_rc.randcnt]);
131}
132
133
134template<UINT32 ALPHA>
135inline void QTIsaac<ALPHA>::randinit(randctx* ctx, bool bUseSeed)
136{
137	UINT32 a,b,c,d,e,f,g,h;
138	int i;
139
140	a = b = c = d = e = f = g = h = GOLDEN_RATIO;
141
142	UINT32* m = (ctx->randmem);
143	UINT32* r = (ctx->randrsl);
144
145	if(!bUseSeed)
146	{
147		ctx->randa = 0;
148		ctx->randb = 0;
149		ctx->randc = 0;
150	}
151
152	// scramble it
153	for(i=0; i < 4; ++i)
154	{
155		shuffle(a,b,c,d,e,f,g,h);
156	}
157
158	if(bUseSeed)
159	{
160		// initialize using the contents of r[] as the seed
161
162		for(i=0; i < N; i+=8)
163		{
164			a+=r[i  ]; b+=r[i+1]; c+=r[i+2]; d+=r[i+3];
165			e+=r[i+4]; f+=r[i+5]; g+=r[i+6]; h+=r[i+7];
166
167			shuffle(a,b,c,d,e,f,g,h);
168
169			m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
170			m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
171		}
172
173		//do a second pass to make all of the seed affect all of m
174
175		for(i=0; i < N; i += 8)
176		{
177			a+=m[i  ]; b+=m[i+1]; c+=m[i+2]; d+=m[i+3];
178			e+=m[i+4]; f+=m[i+5]; g+=m[i+6]; h+=m[i+7];
179
180			shuffle(a,b,c,d,e,f,g,h);
181
182			m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
183			m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
184		}
185	}
186	else
187	{
188		// fill in mm[] with messy stuff
189
190		shuffle(a,b,c,d,e,f,g,h);
191
192		m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
193		m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
194
195	}
196
197	isaac(ctx);         // fill in the first set of results
198	ctx->randcnt = N;   // prepare to use the first set of results
199}
200
201
202template<UINT32 ALPHA>
203inline UINT32 QTIsaac<ALPHA>::ind(UINT32* mm, UINT32 x)
204{
205	return (*(UINT32*)((byte*)(mm) + ((x) & ((N-1)<<2))));
206}
207
208
209template<UINT32 ALPHA>
210inline void QTIsaac<ALPHA>::rngstep(UINT32 mix, UINT32& a, UINT32& b, UINT32*& mm, UINT32*& m, UINT32*& m2, UINT32*& r, UINT32& x, UINT32& y)
211{
212	x = *m;
213	a = (a^(mix)) + *(m2++);
214	*(m++) = y = ind(mm,x) + a + b;
215	*(r++) = b = ind(mm,y>>ALPHA) + x;
216}
217
218
219template<UINT32 ALPHA>
220inline void QTIsaac<ALPHA>::shuffle(UINT32& a, UINT32& b, UINT32& c, UINT32& d, UINT32& e, UINT32& f, UINT32& g, UINT32& h)
221{
222	a^=b<<11; d+=a; b+=c;
223	b^=c>>2;  e+=b; c+=d;
224	c^=d<<8;  f+=c; d+=e;
225	d^=e>>16; g+=d; e+=f;
226	e^=f<<10; h+=e; f+=g;
227	f^=g>>4;  a+=f; g+=h;
228	g^=h<<8;  b+=g; h+=a;
229	h^=a>>9;  c+=h; a+=b;
230}
231
232
233template<UINT32 ALPHA>
234inline void QTIsaac<ALPHA>::isaac(randctx* ctx)
235{
236	UINT32 x,y;
237
238	UINT32* mm = ctx->randmem;
239	UINT32* r  = ctx->randrsl;
240
241	UINT32 a = (ctx->randa);
242	UINT32 b = (ctx->randb + (++ctx->randc));
243
244	UINT32* m    = mm;
245	UINT32* m2   = (m+(N/2));
246	UINT32* mend = m2;
247
248	for(; m<mend; )
249	{
250		rngstep((a<<13), a, b, mm, m, m2, r, x, y);
251		rngstep((a>>6) , a, b, mm, m, m2, r, x, y);
252		rngstep((a<<2) , a, b, mm, m, m2, r, x, y);
253		rngstep((a>>16), a, b, mm, m, m2, r, x, y);
254	}
255
256	m2 = mm;
257
258	for(; m2<mend; )
259	{
260		rngstep((a<<13), a, b, mm, m, m2, r, x, y);
261		rngstep((a>>6) , a, b, mm, m, m2, r, x, y);
262		rngstep((a<<2) , a, b, mm, m, m2, r, x, y);
263		rngstep((a>>16), a, b, mm, m, m2, r, x, y);
264	}
265
266	ctx->randb = b;
267	ctx->randa = a;
268}
269
270
271#endif // __ISAAC_HPP
272
273