Tim's ISAAC module.
[ccan] / ccan / isaac / isaac.c
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.*/
3 #include <float.h>
4 #include <math.h>
5 #include <string.h>
6 #include <ccan/ilog/ilog.h>
7 #include "isaac.h"
8 #if defined(__GNUC_PREREQ)
9 # if __GNUC_PREREQ(4,2)
10 #  pragma GCC diagnostic ignored "-Wparentheses"
11 # endif
12 #endif
13
14
15
16 #define ISAAC_MASK        (0xFFFFFFFFU)
17
18
19
20 static void isaac_update(isaac_ctx *_ctx){
21   uint32_t *m;
22   uint32_t *r;
23   uint32_t  a;
24   uint32_t  b;
25   uint32_t  x;
26   uint32_t  y;
27   int       i;
28   m=_ctx->m;
29   r=_ctx->r;
30   a=_ctx->a;
31   b=_ctx->b+(++_ctx->c)&ISAAC_MASK;
32   for(i=0;i<ISAAC_SZ/2;i++){
33     x=m[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;
37     x=m[++i];
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;
41     x=m[++i];
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;
45     x=m[++i];
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;
49   }
50   for(i=ISAAC_SZ/2;i<ISAAC_SZ;i++){
51     x=m[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;
55     x=m[++i];
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;
59     x=m[++i];
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;
63     x=m[++i];
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;
67   }
68   _ctx->b=b;
69   _ctx->a=a;
70   _ctx->n=ISAAC_SZ;
71 }
72
73 static void isaac_mix(uint32_t _x[8]){
74   static const unsigned char SHIFT[8]={11,2,8,16,10,4,8,9};
75   int i;
76   for(i=0;i<8;i++){
77     _x[i]^=_x[i+1&7]<<SHIFT[i];
78     _x[i+3&7]+=_x[i];
79     _x[i+1&7]+=_x[i+2&7];
80     i++;
81     _x[i]^=_x[i+1&7]>>SHIFT[i];
82     _x[i+3&7]+=_x[i];
83     _x[i+1&7]+=_x[i+2&7];
84   }
85 }
86
87
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);
92 }
93
94 void isaac_reseed(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
95   uint32_t *m;
96   uint32_t *r;
97   uint32_t  x[8];
98   int       i;
99   int       j;
100   m=_ctx->m;
101   r=_ctx->r;
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];
106   }
107   _nseed-=i<<2;
108   if(_nseed>0){
109     uint32_t ri;
110     ri=_seed[i<<2];
111     for(j=1;j<_nseed;j++)ri|=(uint32_t)_seed[i<<2|j]<<(j<<3);
112     r[i++]^=ri;
113   }
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];
118     isaac_mix(x);
119     memcpy(m+i,x,sizeof(x));
120   }
121   for(i=0;i<ISAAC_SZ;i+=8){
122     for(j=0;j<8;j++)x[j]+=m[i+j];
123     isaac_mix(x);
124     memcpy(m+i,x,sizeof(x));
125   }
126   isaac_update(_ctx);
127 }
128
129 uint32_t isaac_next_uint32(isaac_ctx *_ctx){
130   if(!_ctx->n)isaac_update(_ctx);
131   return _ctx->r[--_ctx->n];
132 }
133
134 uint32_t isaac_next_uint(isaac_ctx *_ctx,uint32_t _n){
135   uint32_t r;
136   uint32_t v;
137   uint32_t d;
138   do{
139     r=isaac_next_uint32(_ctx);
140     v=r%_n;
141     d=r-v;
142   }
143   while((d+_n-1&ISAAC_MASK)<d);
144   return v;
145 }
146
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
152            (exclusive).
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){
156   float ret;
157   int   nbits_needed;
158   while(!_bits){
159     if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
160     _base-=32;
161     _bits=isaac_next_uint32(_ctx);
162   }
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);
170 #if FLT_MANT_DIG>32
171   ret=ldexpf((float)_bits,_base);
172 # if FLT_MANT_DIG>65
173   while(32-nbits_needed<0){
174 # else
175   if(32-nbits_needed<0){
176 # endif
177     _base-=32;
178     nbits_needed-=32;
179     ret+=ldexpf((float)isaac_next_uint32(_ctx),_base);
180   }
181   _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
182   ret+=ldexpf((float)_bits,_base-nbits_needed);
183 #else
184   if(nbits_needed>0){
185     _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
186   }
187 # if FLT_MANT_DIG<32
188   else _bits>>=-nbits_needed;
189 # endif
190   ret=ldexpf((float)_bits,_base-nbits_needed);
191 #endif
192   return ret;
193 }
194
195 float isaac_next_float(isaac_ctx *_ctx){
196   return isaac_float_bits(_ctx,0,0);
197 }
198
199 float isaac_next_signed_float(isaac_ctx *_ctx){
200   uint32_t bits;
201   bits=isaac_next_uint32(_ctx);
202   return (1|-((int)bits&1))*isaac_float_bits(_ctx,bits>>1,-31);
203 }
204
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
209            (exclusive).
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){
213   double ret;
214   int    nbits_needed;
215   while(!_bits){
216     if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
217     _base-=32;
218     _bits=isaac_next_uint32(_ctx);
219   }
220   nbits_needed=DBL_MANT_DIG-ILOGNZ_32(_bits);
221 #if DBL_MANT_DIG>32
222   ret=ldexp((double)_bits,_base);
223 # if DBL_MANT_DIG>65
224   while(32-nbits_needed<0){
225 # else
226   if(32-nbits_needed<0){
227 # endif
228     _base-=32;
229     nbits_needed-=32;
230     ret+=ldexp((double)isaac_next_uint32(_ctx),_base);
231   }
232   _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
233   ret+=ldexp((double)_bits,_base-nbits_needed);
234 #else
235   if(nbits_needed>0){
236     _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
237   }
238 # if DBL_MANT_DIG<32
239   else _bits>>=-nbits_needed;
240 # endif
241   ret=ldexp((double)_bits,exp-DBL_MANT_DIG);
242 #endif
243   return ret;
244 }
245
246 double isaac_next_double(isaac_ctx *_ctx){
247   return isaac_double_bits(_ctx,0,0);
248 }
249
250 double isaac_next_signed_double(isaac_ctx *_ctx){
251   uint32_t bits;
252   bits=isaac_next_uint32(_ctx);
253   return (1|-((int)bits&1))*isaac_double_bits(_ctx,bits>>1,-31);
254 }