1 /*Written by Timothy B. Terriberry (tterribe@xiph.org) 1999-2009 public domain.
2 Based on the public domain ISAAC 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 ISAAC64_MASK ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
20 static void isaac64_update(isaac64_ctx *_ctx){
31 b=_ctx->b+(++_ctx->c)&ISAAC64_MASK;
32 for(i=0;i<ISAAC64_SZ/2;i++){
34 a=~(a^a<<21)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
35 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
36 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
38 a=(a^a>>5)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
39 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
40 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
42 a=(a^a<<12)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
43 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
44 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
46 a=(a^a>>33)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
47 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
48 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
50 for(i=ISAAC64_SZ/2;i<ISAAC64_SZ;i++){
52 a=~(a^a<<21)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
53 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
54 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
56 a=(a^a>>5)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
57 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
58 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
60 a=(a^a<<12)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
61 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
62 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
64 a=(a^a>>33)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
65 m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
66 r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
73 static void isaac64_mix(uint64_t _x[8]){
74 static const unsigned char SHIFT[8]={9,9,23,15,14,20,17,14};
78 _x[i+5&7]^=_x[i+7&7]>>SHIFT[i];
82 _x[i+5&7]^=_x[i+7&7]<<SHIFT[i];
88 void isaac64_init(isaac64_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 isaac64_reseed(_ctx,_seed,_nseed);
94 void isaac64_reseed(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed){
102 if(_nseed>ISAAC64_SEED_SZ_MAX)_nseed=ISAAC64_SEED_SZ_MAX;
103 for(i=0;i<_nseed>>3;i++){
104 r[i]^=(uint64_t)_seed[i<<3|7]<<56|(uint64_t)_seed[i<<3|6]<<48|
105 (uint64_t)_seed[i<<3|5]<<40|(uint64_t)_seed[i<<3|4]<<32|
106 (uint64_t)_seed[i<<3|3]<<24|(uint64_t)_seed[i<<3|2]<<16|
107 (uint64_t)_seed[i<<3|1]<<8|_seed[i<<3];
113 for(j=1;j<_nseed;j++)ri|=(uint64_t)_seed[i<<3|j]<<(j<<3);
116 x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=(uint64_t)0x9E3779B97F4A7C13ULL;
117 for(i=0;i<4;i++)isaac64_mix(x);
118 for(i=0;i<ISAAC64_SZ;i+=8){
119 for(j=0;j<8;j++)x[j]+=r[i+j];
121 memcpy(m+i,x,sizeof(x));
123 for(i=0;i<ISAAC64_SZ;i+=8){
124 for(j=0;j<8;j++)x[j]+=m[i+j];
126 memcpy(m+i,x,sizeof(x));
128 isaac64_update(_ctx);
131 uint64_t isaac64_next_uint64(isaac64_ctx *_ctx){
132 if(!_ctx->n)isaac64_update(_ctx);
133 return _ctx->r[--_ctx->n];
136 uint64_t isaac64_next_uint(isaac64_ctx *_ctx,uint64_t _n){
141 r=isaac64_next_uint64(_ctx);
145 while((d+_n-1&ISAAC64_MASK)<d);
149 /*Returns a uniform random float.
150 The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
151 _bits: An initial set of random bits.
152 _base: This should be -(the number of bits in _bits), up to -64.
153 Return: A float uniformly distributed between 0 (inclusive) and 1
155 The average value was measured over 2**32 samples to be
156 0.499991407275206357.*/
157 static float isaac64_float_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
161 if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
163 _bits=isaac64_next_uint64(_ctx);
165 nbits_needed=FLT_MANT_DIG-ILOGNZ_64(_bits);
167 ret=ldexpf((float)_bits,_base);
168 # if FLT_MANT_DIG>129
169 while(64-nbits_needed<0){
171 if(64-nbits_needed<0){
175 ret+=ldexpf((float)isaac64_next_uint64(_ctx),_base);
177 _bits=isaac64_next_uint64(_ctx)>>64-nbits_needed;
178 ret+=ldexpf((float)_bits,_base-nbits_needed);
181 _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>64-nbits_needed;
184 else _bits>>=-nbits_needed;
186 ret=ldexpf((float)_bits,_base-nbits_needed);
191 float isaac64_next_float(isaac64_ctx *_ctx){
192 return isaac64_float_bits(_ctx,0,0);
195 float isaac64_next_signed_float(isaac64_ctx *_ctx){
197 bits=isaac64_next_uint64(_ctx);
198 return (1|-((int)bits&1))*isaac64_float_bits(_ctx,bits>>1,-63);
201 /*Returns a uniform random double.
202 _bits: An initial set of random bits.
203 _base: This should be -(the number of bits in _bits), up to -64.
204 Return: A double uniformly distributed between 0 (inclusive) and 1
206 The average value was measured over 2**32 samples to be
207 0.499990992392019273.*/
208 static double isaac64_double_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
212 if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
214 _bits=isaac64_next_uint64(_ctx);
216 nbits_needed=DBL_MANT_DIG-ILOGNZ_64(_bits);
218 ret=ldexp((double)_bits,_base);
219 # if DBL_MANT_DIG>129
220 while(64-nbits_needed<0){
222 if(64-nbits_needed<0){
226 ret+=ldexp((double)isaac64_next_uint64(_ctx),_base);
228 _bits=isaac64_next_uint64(_ctx)>>64-nbits_needed;
229 ret+=ldexp((double)_bits,_base-nbits_needed);
232 _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>64-nbits_needed;
235 else _bits>>=-nbits_needed;
237 ret=ldexp((double)_bits,_base-nbits_needed);
242 double isaac64_next_double(isaac64_ctx *_ctx){
243 return isaac64_double_bits(_ctx,0,0);
246 double isaac64_next_signed_double(isaac64_ctx *_ctx){
248 bits=isaac64_next_uint64(_ctx);
249 return (1|-((int)bits&1))*isaac64_double_bits(_ctx,bits>>1,-63);