Fix warnings for isaac w/ gcc4.1.
[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
9
10 #define ISAAC_MASK        (0xFFFFFFFFU)
11
12 /* Extract ISAAC_SZ_LOG bits (starting at bit 2). */
13 static inline uint32_t lower_bits(uint32_t x)
14 {
15         return (x & ((ISAAC_SZ-1) << 2)) >> 2;
16 }
17
18 /* Extract next ISAAC_SZ_LOG bits (starting at bit ISAAC_SZ_LOG+2). */
19 static inline uint32_t upper_bits(uint32_t y)
20 {
21         return (y >> (ISAAC_SZ_LOG+2)) & (ISAAC_SZ-1);
22 }
23
24 static void isaac_update(isaac_ctx *_ctx){
25   uint32_t *m;
26   uint32_t *r;
27   uint32_t  a;
28   uint32_t  b;
29   uint32_t  x;
30   uint32_t  y;
31   int       i;
32   m=_ctx->m;
33   r=_ctx->r;
34   a=_ctx->a;
35   b=_ctx->b+(++_ctx->c);
36   for(i=0;i<ISAAC_SZ/2;i++){
37     x=m[i];
38     a=(a^a<<13)+m[i+ISAAC_SZ/2];
39     m[i]=y=m[lower_bits(x)]+a+b;
40     r[i]=b=m[upper_bits(y)]+x;
41     x=m[++i];
42     a=(a^a>>6)+m[i+ISAAC_SZ/2];
43     m[i]=y=m[lower_bits(x)]+a+b;
44     r[i]=b=m[upper_bits(y)]+x;
45     x=m[++i];
46     a=(a^a<<2)+m[i+ISAAC_SZ/2];
47     m[i]=y=m[lower_bits(x)]+a+b;
48     r[i]=b=m[upper_bits(y)]+x;
49     x=m[++i];
50     a=(a^a>>16)+m[i+ISAAC_SZ/2];
51     m[i]=y=m[lower_bits(x)]+a+b;
52     r[i]=b=m[upper_bits(y)]+x;
53   }
54   for(i=ISAAC_SZ/2;i<ISAAC_SZ;i++){
55     x=m[i];
56     a=(a^a<<13)+m[i-ISAAC_SZ/2];
57     m[i]=y=m[lower_bits(x)]+a+b;
58     r[i]=b=m[upper_bits(y)]+x;
59     x=m[++i];
60     a=(a^a>>6)+m[i-ISAAC_SZ/2];
61     m[i]=y=m[lower_bits(x)]+a+b;
62     r[i]=b=m[upper_bits(y)]+x;
63     x=m[++i];
64     a=(a^a<<2)+m[i-ISAAC_SZ/2];
65     m[i]=y=m[lower_bits(x)]+a+b;
66     r[i]=b=m[upper_bits(y)]+x;
67     x=m[++i];
68     a=(a^a>>16)+m[i-ISAAC_SZ/2];
69     m[i]=y=m[lower_bits(x)]+a+b;
70     r[i]=b=m[upper_bits(y)]+x;
71   }
72   _ctx->b=b;
73   _ctx->a=a;
74   _ctx->n=ISAAC_SZ;
75 }
76
77 static void isaac_mix(uint32_t _x[8]){
78   static const unsigned char SHIFT[8]={11,2,8,16,10,4,8,9};
79   int i;
80   for(i=0;i<8;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     i++;
85     _x[i]^=_x[(i+1)&7]>>SHIFT[i];
86     _x[(i+3)&7]+=_x[i];
87     _x[(i+1)&7]+=_x[(i+2)&7];
88   }
89 }
90
91
92 void isaac_init(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
93   _ctx->a=_ctx->b=_ctx->c=0;
94   memset(_ctx->r,0,sizeof(_ctx->r));
95   isaac_reseed(_ctx,_seed,_nseed);
96 }
97
98 void isaac_reseed(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
99   uint32_t *m;
100   uint32_t *r;
101   uint32_t  x[8];
102   int       i;
103   int       j;
104   m=_ctx->m;
105   r=_ctx->r;
106   if(_nseed>ISAAC_SEED_SZ_MAX)_nseed=ISAAC_SEED_SZ_MAX;
107   for(i=0;i<_nseed>>2;i++){
108     r[i]^=(uint32_t)_seed[i<<2|3]<<24|(uint32_t)_seed[i<<2|2]<<16|
109      (uint32_t)_seed[i<<2|1]<<8|_seed[i<<2];
110   }
111   _nseed-=i<<2;
112   if(_nseed>0){
113     uint32_t ri;
114     ri=_seed[i<<2];
115     for(j=1;j<_nseed;j++)ri|=(uint32_t)_seed[i<<2|j]<<(j<<3);
116     r[i++]^=ri;
117   }
118   x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=0x9E3779B9U;
119   for(i=0;i<4;i++)isaac_mix(x);
120   for(i=0;i<ISAAC_SZ;i+=8){
121     for(j=0;j<8;j++)x[j]+=r[i+j];
122     isaac_mix(x);
123     memcpy(m+i,x,sizeof(x));
124   }
125   for(i=0;i<ISAAC_SZ;i+=8){
126     for(j=0;j<8;j++)x[j]+=m[i+j];
127     isaac_mix(x);
128     memcpy(m+i,x,sizeof(x));
129   }
130   isaac_update(_ctx);
131 }
132
133 uint32_t isaac_next_uint32(isaac_ctx *_ctx){
134   if(!_ctx->n)isaac_update(_ctx);
135   return _ctx->r[--_ctx->n];
136 }
137
138 uint32_t isaac_next_uint(isaac_ctx *_ctx,uint32_t _n){
139   uint32_t r;
140   uint32_t v;
141   uint32_t d;
142   do{
143     r=isaac_next_uint32(_ctx);
144     v=r%_n;
145     d=r-v;
146   }
147   while(((d+_n-1)&ISAAC_MASK)<d);
148   return v;
149 }
150
151 /*Returns a uniform random float.
152   The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
153   _bits: An initial set of random bits.
154   _base: This should be -(the number of bits in _bits), up to -32.
155   Return: A float uniformly distributed between 0 (inclusive) and 1
156            (exclusive).
157           The average value was measured over 2**32 samples to be
158            0.50000037448772916.*/
159 static float isaac_float_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
160   float ret;
161   int   nbits_needed;
162   while(!_bits){
163     if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
164     _base-=32;
165     _bits=isaac_next_uint32(_ctx);
166   }
167   /*Note: This could also be determined with frexp(), for a slightly more
168      portable solution, but that takes twice as long, and one has to worry
169      about rounding effects, which can over-estimate the exponent when given
170      FLT_MANT_DIG+1 consecutive one bits.
171     Even the fallback C implementation of ILOGNZ_32() yields an implementation
172      25% faster than the frexp() method.*/
173   nbits_needed=FLT_MANT_DIG-ILOGNZ_32(_bits);
174 #if FLT_MANT_DIG>32
175   ret=ldexpf((float)_bits,_base);
176 # if FLT_MANT_DIG>65
177   while(32-nbits_needed<0){
178 # else
179   if(32-nbits_needed<0){
180 # endif
181     _base-=32;
182     nbits_needed-=32;
183     ret+=ldexpf((float)isaac_next_uint32(_ctx),_base);
184   }
185   _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
186   ret+=ldexpf((float)_bits,_base-nbits_needed);
187 #else
188   if(nbits_needed>0){
189     _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>(32-nbits_needed);
190   }
191 # if FLT_MANT_DIG<32
192   else _bits>>=-nbits_needed;
193 # endif
194   ret=ldexpf((float)_bits,_base-nbits_needed);
195 #endif
196   return ret;
197 }
198
199 float isaac_next_float(isaac_ctx *_ctx){
200   return isaac_float_bits(_ctx,0,0);
201 }
202
203 float isaac_next_signed_float(isaac_ctx *_ctx){
204   uint32_t bits;
205   bits=isaac_next_uint32(_ctx);
206   return (1|-((int)bits&1))*isaac_float_bits(_ctx,bits>>1,-31);
207 }
208
209 /*Returns a uniform random double.
210   _bits: An initial set of random bits.
211   _base: This should be -(the number of bits in _bits), up to -32.
212   Return: A double uniformly distributed between 0 (inclusive) and 1
213            (exclusive).
214           The average value was measured over 2**32 samples to be
215            0.500006289408060911*/
216 static double isaac_double_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
217   double ret;
218   int    nbits_needed;
219   while(!_bits){
220     if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
221     _base-=32;
222     _bits=isaac_next_uint32(_ctx);
223   }
224   nbits_needed=DBL_MANT_DIG-ILOGNZ_32(_bits);
225 #if DBL_MANT_DIG>32
226   ret=ldexp((double)_bits,_base);
227 # if DBL_MANT_DIG>65
228   while(32-nbits_needed<0){
229 # else
230   if(32-nbits_needed<0){
231 # endif
232     _base-=32;
233     nbits_needed-=32;
234     ret+=ldexp((double)isaac_next_uint32(_ctx),_base);
235   }
236   _bits=isaac_next_uint32(_ctx)>>(32-nbits_needed);
237   ret+=ldexp((double)_bits,_base-nbits_needed);
238 #else
239   if(nbits_needed>0){
240     _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
241   }
242 # if DBL_MANT_DIG<32
243   else _bits>>=-nbits_needed;
244 # endif
245   ret=ldexp((double)_bits,exp-DBL_MANT_DIG);
246 #endif
247   return ret;
248 }
249
250 double isaac_next_double(isaac_ctx *_ctx){
251   return isaac_double_bits(_ctx,0,0);
252 }
253
254 double isaac_next_signed_double(isaac_ctx *_ctx){
255   uint32_t bits;
256   bits=isaac_next_uint32(_ctx);
257   return (1|-((int)bits&1))*isaac_double_bits(_ctx,bits>>1,-31);
258 }