1 /*Written by Timothy B. Terriberry (tterribe@xiph.org) 1999-2009.
2 CC0 (Public domain) - see LICENSE file for details
3 Based on the public domain implementation by Robert J. Jenkins Jr.*/
7 #include <ccan/ilog/ilog.h>
11 #define ISAAC_MASK (0xFFFFFFFFU)
13 /* Extract ISAAC_SZ_LOG bits (starting at bit 2). */
14 static inline uint32_t lower_bits(uint32_t x)
16 return (x & ((ISAAC_SZ-1) << 2)) >> 2;
19 /* Extract next ISAAC_SZ_LOG bits (starting at bit ISAAC_SZ_LOG+2). */
20 static inline uint32_t upper_bits(uint32_t y)
22 return (y >> (ISAAC_SZ_LOG+2)) & (ISAAC_SZ-1);
25 static void isaac_update(isaac_ctx *_ctx){
36 b=_ctx->b+(++_ctx->c);
37 for(i=0;i<ISAAC_SZ/2;i++){
39 a=(a^a<<13)+m[i+ISAAC_SZ/2];
40 m[i]=y=m[lower_bits(x)]+a+b;
41 r[i]=b=m[upper_bits(y)]+x;
43 a=(a^a>>6)+m[i+ISAAC_SZ/2];
44 m[i]=y=m[lower_bits(x)]+a+b;
45 r[i]=b=m[upper_bits(y)]+x;
47 a=(a^a<<2)+m[i+ISAAC_SZ/2];
48 m[i]=y=m[lower_bits(x)]+a+b;
49 r[i]=b=m[upper_bits(y)]+x;
51 a=(a^a>>16)+m[i+ISAAC_SZ/2];
52 m[i]=y=m[lower_bits(x)]+a+b;
53 r[i]=b=m[upper_bits(y)]+x;
55 for(i=ISAAC_SZ/2;i<ISAAC_SZ;i++){
57 a=(a^a<<13)+m[i-ISAAC_SZ/2];
58 m[i]=y=m[lower_bits(x)]+a+b;
59 r[i]=b=m[upper_bits(y)]+x;
61 a=(a^a>>6)+m[i-ISAAC_SZ/2];
62 m[i]=y=m[lower_bits(x)]+a+b;
63 r[i]=b=m[upper_bits(y)]+x;
65 a=(a^a<<2)+m[i-ISAAC_SZ/2];
66 m[i]=y=m[lower_bits(x)]+a+b;
67 r[i]=b=m[upper_bits(y)]+x;
69 a=(a^a>>16)+m[i-ISAAC_SZ/2];
70 m[i]=y=m[lower_bits(x)]+a+b;
71 r[i]=b=m[upper_bits(y)]+x;
78 static void isaac_mix(uint32_t _x[8]){
79 static const unsigned char SHIFT[8]={11,2,8,16,10,4,8,9};
82 _x[i]^=_x[(i+1)&7]<<SHIFT[i];
84 _x[(i+1)&7]+=_x[(i+2)&7];
86 _x[i]^=_x[(i+1)&7]>>SHIFT[i];
88 _x[(i+1)&7]+=_x[(i+2)&7];
93 void isaac_init(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
94 _ctx->a=_ctx->b=_ctx->c=0;
95 memset(_ctx->r,0,sizeof(_ctx->r));
96 isaac_reseed(_ctx,_seed,_nseed);
99 void isaac_reseed(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
107 if(_nseed>ISAAC_SEED_SZ_MAX)_nseed=ISAAC_SEED_SZ_MAX;
108 for(i=0;i<_nseed>>2;i++){
109 r[i]^=(uint32_t)_seed[i<<2|3]<<24|(uint32_t)_seed[i<<2|2]<<16|
110 (uint32_t)_seed[i<<2|1]<<8|_seed[i<<2];
116 for(j=1;j<_nseed;j++)ri|=(uint32_t)_seed[i<<2|j]<<(j<<3);
119 x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=0x9E3779B9U;
120 for(i=0;i<4;i++)isaac_mix(x);
121 for(i=0;i<ISAAC_SZ;i+=8){
122 for(j=0;j<8;j++)x[j]+=r[i+j];
124 memcpy(m+i,x,sizeof(x));
126 for(i=0;i<ISAAC_SZ;i+=8){
127 for(j=0;j<8;j++)x[j]+=m[i+j];
129 memcpy(m+i,x,sizeof(x));
134 uint32_t isaac_next_uint32(isaac_ctx *_ctx){
135 if(!_ctx->n)isaac_update(_ctx);
136 return _ctx->r[--_ctx->n];
139 uint32_t isaac_next_uint(isaac_ctx *_ctx,uint32_t _n){
144 r=isaac_next_uint32(_ctx);
148 while(((d+_n-1)&ISAAC_MASK)<d);
152 /*Returns a uniform random float.
153 The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
154 _bits: An initial set of random bits.
155 _base: This should be -(the number of bits in _bits), up to -32.
156 Return: A float uniformly distributed between 0 (inclusive) and 1
158 The average value was measured over 2**32 samples to be
159 0.50000037448772916.*/
160 static float isaac_float_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
164 if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
166 _bits=isaac_next_uint32(_ctx);
168 /*Note: This could also be determined with frexp(), for a slightly more
169 portable solution, but that takes twice as long, and one has to worry
170 about rounding effects, which can over-estimate the exponent when given
171 FLT_MANT_DIG+1 consecutive one bits.
172 Even the fallback C implementation of ILOGNZ_32() yields an implementation
173 25% faster than the frexp() method.*/
174 nbits_needed=FLT_MANT_DIG-ilog32_nz(_bits);
176 ret=ldexpf((float)_bits,_base);
178 while(32-nbits_needed<0){
180 if(32-nbits_needed<0){
184 ret+=ldexpf((float)isaac_next_uint32(_ctx),_base);
186 _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
187 ret+=ldexpf((float)_bits,_base-nbits_needed);
190 _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>(32-nbits_needed);
193 else _bits>>=-nbits_needed;
195 ret=ldexpf((float)_bits,_base-nbits_needed);
200 float isaac_next_float(isaac_ctx *_ctx){
201 return isaac_float_bits(_ctx,0,0);
204 float isaac_next_signed_float(isaac_ctx *_ctx){
206 bits=isaac_next_uint32(_ctx);
207 return (1|-((int)bits&1))*isaac_float_bits(_ctx,bits>>1,-31);
210 /*Returns a uniform random double.
211 _bits: An initial set of random bits.
212 _base: This should be -(the number of bits in _bits), up to -32.
213 Return: A double uniformly distributed between 0 (inclusive) and 1
215 The average value was measured over 2**32 samples to be
216 0.500006289408060911*/
217 static double isaac_double_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
221 if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
223 _bits=isaac_next_uint32(_ctx);
225 nbits_needed=DBL_MANT_DIG-ilog32_nz(_bits);
227 ret=ldexp((double)_bits,_base);
229 while(32-nbits_needed<0){
231 if(32-nbits_needed<0){
235 ret+=ldexp((double)isaac_next_uint32(_ctx),_base);
237 _bits=isaac_next_uint32(_ctx)>>(32-nbits_needed);
238 ret+=ldexp((double)_bits,_base-nbits_needed);
241 _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
244 else _bits>>=-nbits_needed;
246 ret=ldexp((double)_bits,exp-DBL_MANT_DIG);
251 double isaac_next_double(isaac_ctx *_ctx){
252 return isaac_double_bits(_ctx,0,0);
255 double isaac_next_signed_double(isaac_ctx *_ctx){
257 bits=isaac_next_uint32(_ctx);
258 return (1|-((int)bits&1))*isaac_double_bits(_ctx,bits>>1,-31);