1 /* rand.c: The ISAAC Pseudo-Random Number Generator
2 *
3 * This is the original ISAAC reference implementation, written by Bob Jenkins
4 * and released into the public domain. The original code by Bob Jenkins was
5 * retrieved from: http://burtleburtle.net/bob/rand/isaacafa.html
6 *
7 * Original filename was rand.c and carried this changelog:
8 * 960327: Creation (addition of randinit, really)
9 * 970719: use context, not global variables, for internal state
10 * 980324: make a portable version
11 * 010626: Note this is public domain
12 *
13 * Jonathan Yu <jawnsy@cpan.org> made some mostly cosmetic changes and
14 * prepared the file for life as a CPAN XS module.
15 *
16 * This package and its contents are released by the author into the
17 * Public Domain, to the full extent permissible by law. For additional
18 * information, please see the included `LICENSE' file.
19 */
20
21 #include "standard.h"
22 #include "rand.h"
23
24 #ifdef USE_PORTABLE
25 #define cut(a) ((a) & 0xffffffff) /* Cut the integer down to 32bits */
26 #else
27 #define cut(a) (a) /* A no-op */
28 #endif
29
30 #define ind(mm,x) ((mm)[(x>>2)&(RANDSIZ-1)])
31 /* the call to cut() is a macro defined in standard.h */
32 #define rngstep(mix,a,b,mm,m,m2,r,x) \
33 { \
34 x = *m; \
35 a = cut((a^(mix)) + *(m2++)); \
36 *(m++) = y = cut(ind(mm,x) + a + b); \
37 *(r++) = b = cut(ind(mm,y>>RANDSIZL) + x); \
38 }
39 #define mix(a,b,c,d,e,f,g,h) \
40 { \
41 a^=b<<11; d+=a; b+=c; \
42 b^=c>>2; e+=b; c+=d; \
43 c^=d<<8; f+=c; d+=e; \
44 d^=e>>16; g+=d; e+=f; \
45 e^=f<<10; h+=e; f+=g; \
46 f^=g>>4; a+=f; g+=h; \
47 g^=h<<8; b+=g; h+=a; \
48 h^=a>>9; c+=h; a+=b; \
49 }
50 #define shuffle(a, b, mm, m, m2, r, x) \
51 { \
52 rngstep(a<<13, a, b, mm, m, m2, r, x); \
53 rngstep(a>>6 , a, b, mm, m, m2, r, x); \
54 rngstep(a<<2 , a, b, mm, m, m2, r, x); \
55 rngstep(a>>16, a, b, mm, m, m2, r, x); \
56 }
57
isaac(randctx * ctx)58 static void isaac(randctx *ctx)
59 {
60 /* Keep these in CPU registers if possible, for speed */
61 register ub4 a, b, x, y;
62 register ub4 *m, *mm, *m2, *r, *mend;
63
64 mm = ctx->randmem;
65 r = ctx->randrsl;
66 a = ctx->randa;
67 b = cut(ctx->randb + (++ctx->randc));
68
69 m = mm;
70 mend = m2 = m + (RANDSIZ / 2);
71 while (m < mend) {
72 shuffle(a, b, mm, m, m2, r, x);
73 }
74
75 m2 = mm;
76 while (m2 < mend) {
77 shuffle(a, b, mm, m, m2, r, x);
78 }
79
80 ctx->randb = b;
81 ctx->randa = a;
82 }
83
84 /* using randrsl[0..RANDSIZ-1] as the seed */
randinit(randctx * ctx)85 void randinit(randctx *ctx)
86 {
87 ub4 a, b, c, d, e, f, g, h;
88 ub4 *m, *r;
89 int i; /* for loop incrementing variable */
90
91 m = ctx->randmem;
92 r = ctx->randrsl;
93
94 ctx->randa = ctx->randb = ctx->randc = (ub4)0;
95
96 /* Initialize a to h with the golden ratio */
97 a=b=c=d=e=f=g=h = 0x9e3779b9;
98
99 /* scramble it */
100 for (i = 0; i < 4; i++) {
101 mix(a,b,c,d,e,f,g,h);
102 }
103
104 /* initialize using the contents of r[] as the seed */
105 for (i = 0; i < RANDSIZ; i += 8) {
106 a += r[i ];
107 b += r[i+1];
108 c += r[i+2];
109 d += r[i+3];
110 e += r[i+4];
111 f += r[i+5];
112 g += r[i+6];
113 h += r[i+7];
114
115 mix(a,b,c,d,e,f,g,h);
116
117 m[i ] = a;
118 m[i+1] = b;
119 m[i+2] = c;
120 m[i+3] = d;
121 m[i+4] = e;
122 m[i+5] = f;
123 m[i+6] = g;
124 m[i+7] = h;
125 }
126
127 /* do a second pass to make all of the seed affect all of m */
128 for (i = 0; i < RANDSIZ; i += 8) {
129 a += m[i ];
130 b += m[i+1];
131 c += m[i+2];
132 d += m[i+3];
133 e += m[i+4];
134 f += m[i+5];
135 g += m[i+6];
136 h += m[i+7];
137
138 mix(a,b,c,d,e,f,g,h);
139
140 m[i ] = a;
141 m[i+1] = b;
142 m[i+2] = c;
143 m[i+3] = d;
144 m[i+4] = e;
145 m[i+5] = f;
146 m[i+6] = g;
147 m[i+7] = h;
148 }
149
150 isaac(ctx); /* fill in the first set of results */
151 ctx->randcnt = RANDSIZ; /* prepare to use the first set of results */
152 }
153
154 /* This function was added by Jonathan Yu to return the next integer (taking
155 * the code out of a macro and putting it into a function instead
156 */
randInt(randctx * ctx)157 ub4 randInt(randctx *ctx)
158 {
159 /* If we run out of numbers, reset the sequence */
160 if (ctx->randcnt-- == 0) {
161 isaac(ctx);
162 ctx->randcnt = RANDSIZ - 1;
163 }
164 return ctx->randrsl[ctx->randcnt];
165 }
166