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 ISAAC implementation by Robert J. Jenkins Jr.*/
7 #include <ccan/ilog/ilog.h>
11 #define ISAAC64_MASK ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
13 /* Extract ISAAC64_SZ_LOG bits (starting at bit 3). */
14 static inline uint32_t lower_bits(uint64_t x)
16 return (x & ((ISAAC64_SZ-1) << 3)) >>3;
19 /* Extract next ISAAC64_SZ_LOG bits (starting at bit ISAAC64_SZ_LOG+2). */
20 static inline uint32_t upper_bits(uint32_t y)
22 return (y >> (ISAAC64_SZ_LOG+3)) & (ISAAC64_SZ-1);
25 static void isaac64_update(isaac64_ctx *_ctx){
36 b=_ctx->b+(++_ctx->c);
37 for(i=0;i<ISAAC64_SZ/2;i++){
39 a=~(a^a<<21)+m[i+ISAAC64_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>>5)+m[i+ISAAC64_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<<12)+m[i+ISAAC64_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>>33)+m[i+ISAAC64_SZ/2];
52 m[i]=y=m[lower_bits(x)]+a+b;
53 r[i]=b=m[upper_bits(y)]+x;
55 for(i=ISAAC64_SZ/2;i<ISAAC64_SZ;i++){
57 a=~(a^a<<21)+m[i-ISAAC64_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>>5)+m[i-ISAAC64_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<<12)+m[i-ISAAC64_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>>33)+m[i-ISAAC64_SZ/2];
70 m[i]=y=m[lower_bits(x)]+a+b;
71 r[i]=b=m[upper_bits(y)]+x;
78 static void isaac64_mix(uint64_t _x[8]){
79 static const unsigned char SHIFT[8]={9,9,23,15,14,20,17,14};
83 _x[(i+5)&7]^=_x[(i+7)&7]>>SHIFT[i];
87 _x[(i+5)&7]^=_x[(i+7)&7]<<SHIFT[i];
93 void isaac64_init(isaac64_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 isaac64_reseed(_ctx,_seed,_nseed);
99 void isaac64_reseed(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed){
107 if(_nseed>ISAAC64_SEED_SZ_MAX)_nseed=ISAAC64_SEED_SZ_MAX;
108 for(i=0;i<_nseed>>3;i++){
109 r[i]^=(uint64_t)_seed[i<<3|7]<<56|(uint64_t)_seed[i<<3|6]<<48|
110 (uint64_t)_seed[i<<3|5]<<40|(uint64_t)_seed[i<<3|4]<<32|
111 (uint64_t)_seed[i<<3|3]<<24|(uint64_t)_seed[i<<3|2]<<16|
112 (uint64_t)_seed[i<<3|1]<<8|_seed[i<<3];
118 for(j=1;j<_nseed;j++)ri|=(uint64_t)_seed[i<<3|j]<<(j<<3);
121 x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=(uint64_t)0x9E3779B97F4A7C13ULL;
122 for(i=0;i<4;i++)isaac64_mix(x);
123 for(i=0;i<ISAAC64_SZ;i+=8){
124 for(j=0;j<8;j++)x[j]+=r[i+j];
126 memcpy(m+i,x,sizeof(x));
128 for(i=0;i<ISAAC64_SZ;i+=8){
129 for(j=0;j<8;j++)x[j]+=m[i+j];
131 memcpy(m+i,x,sizeof(x));
133 isaac64_update(_ctx);
136 uint64_t isaac64_next_uint64(isaac64_ctx *_ctx){
137 if(!_ctx->n)isaac64_update(_ctx);
138 return _ctx->r[--_ctx->n];
141 uint64_t isaac64_next_uint(isaac64_ctx *_ctx,uint64_t _n){
146 r=isaac64_next_uint64(_ctx);
150 while(((d+_n-1)&ISAAC64_MASK)<d);
154 /*Returns a uniform random float.
155 The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
156 _bits: An initial set of random bits.
157 _base: This should be -(the number of bits in _bits), up to -64.
158 Return: A float uniformly distributed between 0 (inclusive) and 1
160 The average value was measured over 2**32 samples to be
161 0.499991407275206357.*/
162 static float isaac64_float_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
166 if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
168 _bits=isaac64_next_uint64(_ctx);
170 nbits_needed=FLT_MANT_DIG-ilog64_nz(_bits);
172 ret=ldexpf((float)_bits,_base);
173 # if FLT_MANT_DIG>129
174 while(64-nbits_needed<0){
176 if(64-nbits_needed<0){
180 ret+=ldexpf((float)isaac64_next_uint64(_ctx),_base);
182 _bits=isaac64_next_uint64(_ctx)>>(64-nbits_needed);
183 ret+=ldexpf((float)_bits,_base-nbits_needed);
186 _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>(64-nbits_needed);
189 else _bits>>=-nbits_needed;
191 ret=ldexpf((float)_bits,_base-nbits_needed);
196 float isaac64_next_float(isaac64_ctx *_ctx){
197 return isaac64_float_bits(_ctx,0,0);
200 float isaac64_next_signed_float(isaac64_ctx *_ctx){
202 bits=isaac64_next_uint64(_ctx);
203 return (1|-((int)bits&1))*isaac64_float_bits(_ctx,bits>>1,-63);
206 /*Returns a uniform random double.
207 _bits: An initial set of random bits.
208 _base: This should be -(the number of bits in _bits), up to -64.
209 Return: A double uniformly distributed between 0 (inclusive) and 1
211 The average value was measured over 2**32 samples to be
212 0.499990992392019273.*/
213 static double isaac64_double_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
217 if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
219 _bits=isaac64_next_uint64(_ctx);
221 nbits_needed=DBL_MANT_DIG-ilog64_nz(_bits);
223 ret=ldexp((double)_bits,_base);
224 # if DBL_MANT_DIG>129
225 while(64-nbits_needed<0){
227 if(64-nbits_needed<0){
231 ret+=ldexp((double)isaac64_next_uint64(_ctx),_base);
233 _bits=isaac64_next_uint64(_ctx)>>(64-nbits_needed);
234 ret+=ldexp((double)_bits,_base-nbits_needed);
237 _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>(64-nbits_needed);
240 else _bits>>=-nbits_needed;
242 ret=ldexp((double)_bits,_base-nbits_needed);
247 double isaac64_next_double(isaac64_ctx *_ctx){
248 return isaac64_double_bits(_ctx,0,0);
251 double isaac64_next_signed_double(isaac64_ctx *_ctx){
253 bits=isaac64_next_uint64(_ctx);
254 return (1|-((int)bits&1))*isaac64_double_bits(_ctx,bits>>1,-63);