json_out: make json_out_finished finish buffer.
[ccan] / ccan / isaac / isaac64.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 ISAAC 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 "isaac64.h"
9
10
11 #define ISAAC64_MASK ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
12
13 /* Extract ISAAC64_SZ_LOG bits (starting at bit 3). */
14 static inline uint32_t lower_bits(uint64_t x)
15 {
16         return (x & ((ISAAC64_SZ-1) << 3)) >>3;
17 }
18
19 /* Extract next ISAAC64_SZ_LOG bits (starting at bit ISAAC64_SZ_LOG+2). */
20 static inline uint32_t upper_bits(uint32_t y)
21 {
22         return (y >> (ISAAC64_SZ_LOG+3)) & (ISAAC64_SZ-1);
23 }
24
25 static void isaac64_update(isaac64_ctx *_ctx){
26   uint64_t *m;
27   uint64_t *r;
28   uint64_t  a;
29   uint64_t  b;
30   uint64_t  x;
31   uint64_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<ISAAC64_SZ/2;i++){
38     x=m[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;
42     x=m[++i];
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;
46     x=m[++i];
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;
50     x=m[++i];
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;
54   }
55   for(i=ISAAC64_SZ/2;i<ISAAC64_SZ;i++){
56     x=m[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;
60     x=m[++i];
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;
64     x=m[++i];
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;
68     x=m[++i];
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;
72   }
73   _ctx->b=b;
74   _ctx->a=a;
75   _ctx->n=ISAAC64_SZ;
76 }
77
78 static void isaac64_mix(uint64_t _x[8]){
79   static const unsigned char SHIFT[8]={9,9,23,15,14,20,17,14};
80   int i;
81   for(i=0;i<8;i++){
82     _x[i]-=_x[(i+4)&7];
83     _x[(i+5)&7]^=_x[(i+7)&7]>>SHIFT[i];
84     _x[(i+7)&7]+=_x[i];
85     i++;
86     _x[i]-=_x[(i+4)&7];
87     _x[(i+5)&7]^=_x[(i+7)&7]<<SHIFT[i];
88     _x[(i+7)&7]+=_x[i];
89   }
90 }
91
92
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);
97 }
98
99 void isaac64_reseed(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed){
100   uint64_t *m;
101   uint64_t *r;
102   uint64_t  x[8];
103   int       i;
104   int       j;
105   m=_ctx->m;
106   r=_ctx->r;
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];
113   }
114   _nseed-=i<<3;
115   if(_nseed>0){
116     uint64_t ri;
117     ri=_seed[i<<3];
118     for(j=1;j<_nseed;j++)ri|=(uint64_t)_seed[i<<3|j]<<(j<<3);
119     r[i++]^=ri;
120   }
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];
125     isaac64_mix(x);
126     memcpy(m+i,x,sizeof(x));
127   }
128   for(i=0;i<ISAAC64_SZ;i+=8){
129     for(j=0;j<8;j++)x[j]+=m[i+j];
130     isaac64_mix(x);
131     memcpy(m+i,x,sizeof(x));
132   }
133   isaac64_update(_ctx);
134 }
135
136 uint64_t isaac64_next_uint64(isaac64_ctx *_ctx){
137   if(!_ctx->n)isaac64_update(_ctx);
138   return _ctx->r[--_ctx->n];
139 }
140
141 uint64_t isaac64_next_uint(isaac64_ctx *_ctx,uint64_t _n){
142   uint64_t r;
143   uint64_t v;
144   uint64_t d;
145   do{
146     r=isaac64_next_uint64(_ctx);
147     v=r%_n;
148     d=r-v;
149   }
150   while(((d+_n-1)&ISAAC64_MASK)<d);
151   return v;
152 }
153
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
159            (exclusive).
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){
163   float ret;
164   int   nbits_needed;
165   while(!_bits){
166     if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
167     _base-=64;
168     _bits=isaac64_next_uint64(_ctx);
169   }
170   nbits_needed=FLT_MANT_DIG-ilog64_nz(_bits);
171 #if FLT_MANT_DIG>64
172   ret=ldexpf((float)_bits,_base);
173 # if FLT_MANT_DIG>129
174   while(64-nbits_needed<0){
175 # else
176   if(64-nbits_needed<0){
177 # endif
178     _base-=64;
179     nbits_needed-=64;
180     ret+=ldexpf((float)isaac64_next_uint64(_ctx),_base);
181   }
182   _bits=isaac64_next_uint64(_ctx)>>(64-nbits_needed);
183   ret+=ldexpf((float)_bits,_base-nbits_needed);
184 #else
185   if(nbits_needed>0){
186     _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>(64-nbits_needed);
187   }
188 # if FLT_MANT_DIG<64
189   else _bits>>=-nbits_needed;
190 # endif
191   ret=ldexpf((float)_bits,_base-nbits_needed);
192 #endif
193   return ret;
194 }
195
196 float isaac64_next_float(isaac64_ctx *_ctx){
197   return isaac64_float_bits(_ctx,0,0);
198 }
199
200 float isaac64_next_signed_float(isaac64_ctx *_ctx){
201   uint64_t bits;
202   bits=isaac64_next_uint64(_ctx);
203   return (1|-((int)bits&1))*isaac64_float_bits(_ctx,bits>>1,-63);
204 }
205
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
210            (exclusive).
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){
214   double ret;
215   int    nbits_needed;
216   while(!_bits){
217     if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
218     _base-=64;
219     _bits=isaac64_next_uint64(_ctx);
220   }
221   nbits_needed=DBL_MANT_DIG-ilog64_nz(_bits);
222 #if DBL_MANT_DIG>64
223   ret=ldexp((double)_bits,_base);
224 # if DBL_MANT_DIG>129
225   while(64-nbits_needed<0){
226 # else
227   if(64-nbits_needed<0){
228 # endif
229     _base-=64;
230     nbits_needed-=64;
231     ret+=ldexp((double)isaac64_next_uint64(_ctx),_base);
232   }
233   _bits=isaac64_next_uint64(_ctx)>>(64-nbits_needed);
234   ret+=ldexp((double)_bits,_base-nbits_needed);
235 #else
236   if(nbits_needed>0){
237     _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>(64-nbits_needed);
238   }
239 # if DBL_MANT_DIG<64
240   else _bits>>=-nbits_needed;
241 # endif
242   ret=ldexp((double)_bits,_base-nbits_needed);
243 #endif
244   return ret;
245 }
246
247 double isaac64_next_double(isaac64_ctx *_ctx){
248   return isaac64_double_bits(_ctx,0,0);
249 }
250
251 double isaac64_next_signed_double(isaac64_ctx *_ctx){
252   uint64_t bits;
253   bits=isaac64_next_uint64(_ctx);
254   return (1|-((int)bits&1))*isaac64_double_bits(_ctx,bits>>1,-63);
255 }