1 /*Written by Timothy B. Terriberry (tterribe@xiph.org) 1999-2009 public domain.
2 Based on the public domain implementation by Robert J. Jenkins Jr.*/
6 #include <ccan/ilog/ilog.h>
8 #if defined(__GNUC_PREREQ)
9 # if __GNUC_PREREQ(4,2)
10 # pragma GCC diagnostic ignored "-Wparentheses"
16 #define ISAAC_MASK (0xFFFFFFFFU)
20 static void isaac_update(isaac_ctx *_ctx){
31 b=_ctx->b+(++_ctx->c)&ISAAC_MASK;
32 for(i=0;i<ISAAC_SZ/2;i++){
34 a=(a^a<<13)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
35 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
36 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
38 a=(a^a>>6)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
39 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
40 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
42 a=(a^a<<2)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
43 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
44 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
46 a=(a^a>>16)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
47 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
48 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
50 for(i=ISAAC_SZ/2;i<ISAAC_SZ;i++){
52 a=(a^a<<13)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
53 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
54 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
56 a=(a^a>>6)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
57 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
58 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
60 a=(a^a<<2)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
61 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
62 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
64 a=(a^a>>16)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
65 m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
66 r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
73 static void isaac_mix(uint32_t _x[8]){
74 static const unsigned char SHIFT[8]={11,2,8,16,10,4,8,9};
77 _x[i]^=_x[i+1&7]<<SHIFT[i];
81 _x[i]^=_x[i+1&7]>>SHIFT[i];
88 void isaac_init(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
89 _ctx->a=_ctx->b=_ctx->c=0;
90 memset(_ctx->r,0,sizeof(_ctx->r));
91 isaac_reseed(_ctx,_seed,_nseed);
94 void isaac_reseed(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
102 if(_nseed>ISAAC_SEED_SZ_MAX)_nseed=ISAAC_SEED_SZ_MAX;
103 for(i=0;i<_nseed>>2;i++){
104 r[i]^=(uint32_t)_seed[i<<2|3]<<24|(uint32_t)_seed[i<<2|2]<<16|
105 (uint32_t)_seed[i<<2|1]<<8|_seed[i<<2];
111 for(j=1;j<_nseed;j++)ri|=(uint32_t)_seed[i<<2|j]<<(j<<3);
114 x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=0x9E3779B9U;
115 for(i=0;i<4;i++)isaac_mix(x);
116 for(i=0;i<ISAAC_SZ;i+=8){
117 for(j=0;j<8;j++)x[j]+=r[i+j];
119 memcpy(m+i,x,sizeof(x));
121 for(i=0;i<ISAAC_SZ;i+=8){
122 for(j=0;j<8;j++)x[j]+=m[i+j];
124 memcpy(m+i,x,sizeof(x));
129 uint32_t isaac_next_uint32(isaac_ctx *_ctx){
130 if(!_ctx->n)isaac_update(_ctx);
131 return _ctx->r[--_ctx->n];
134 uint32_t isaac_next_uint(isaac_ctx *_ctx,uint32_t _n){
139 r=isaac_next_uint32(_ctx);
143 while((d+_n-1&ISAAC_MASK)<d);
147 /*Returns a uniform random float.
148 The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
149 _bits: An initial set of random bits.
150 _base: This should be -(the number of bits in _bits), up to -32.
151 Return: A float uniformly distributed between 0 (inclusive) and 1
153 The average value was measured over 2**32 samples to be
154 0.50000037448772916.*/
155 static float isaac_float_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
159 if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
161 _bits=isaac_next_uint32(_ctx);
163 /*Note: This could also be determined with frexp(), for a slightly more
164 portable solution, but that takes twice as long, and one has to worry
165 about rounding effects, which can over-estimate the exponent when given
166 FLT_MANT_DIG+1 consecutive one bits.
167 Even the fallback C implementation of ILOGNZ_32() yields an implementation
168 25% faster than the frexp() method.*/
169 nbits_needed=FLT_MANT_DIG-ILOGNZ_32(_bits);
171 ret=ldexpf((float)_bits,_base);
173 while(32-nbits_needed<0){
175 if(32-nbits_needed<0){
179 ret+=ldexpf((float)isaac_next_uint32(_ctx),_base);
181 _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
182 ret+=ldexpf((float)_bits,_base-nbits_needed);
185 _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
188 else _bits>>=-nbits_needed;
190 ret=ldexpf((float)_bits,_base-nbits_needed);
195 float isaac_next_float(isaac_ctx *_ctx){
196 return isaac_float_bits(_ctx,0,0);
199 float isaac_next_signed_float(isaac_ctx *_ctx){
201 bits=isaac_next_uint32(_ctx);
202 return (1|-((int)bits&1))*isaac_float_bits(_ctx,bits>>1,-31);
205 /*Returns a uniform random double.
206 _bits: An initial set of random bits.
207 _base: This should be -(the number of bits in _bits), up to -32.
208 Return: A double uniformly distributed between 0 (inclusive) and 1
210 The average value was measured over 2**32 samples to be
211 0.500006289408060911*/
212 static double isaac_double_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
216 if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
218 _bits=isaac_next_uint32(_ctx);
220 nbits_needed=DBL_MANT_DIG-ILOGNZ_32(_bits);
222 ret=ldexp((double)_bits,_base);
224 while(32-nbits_needed<0){
226 if(32-nbits_needed<0){
230 ret+=ldexp((double)isaac_next_uint32(_ctx),_base);
232 _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
233 ret+=ldexp((double)_bits,_base-nbits_needed);
236 _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
239 else _bits>>=-nbits_needed;
241 ret=ldexp((double)_bits,exp-DBL_MANT_DIG);
246 double isaac_next_double(isaac_ctx *_ctx){
247 return isaac_double_bits(_ctx,0,0);
250 double isaac_next_signed_double(isaac_ctx *_ctx){
252 bits=isaac_next_uint32(_ctx);
253 return (1|-((int)bits&1))*isaac_double_bits(_ctx,bits>>1,-31);