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