Use -O not -O3: reduces ccan/tdb test time from 24 to 18 seconds.
[ccan] / ccan / isaac / isaac64.c
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.*/
3 #include <float.h>
4 #include <math.h>
5 #include <string.h>
6 #include <ccan/ilog/ilog.h>
7 #include "isaac64.h"
8
9
10 #define ISAAC64_MASK ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
11
12 /* Extract ISAAC64_SZ_LOG bits (starting at bit 3). */
13 static inline uint32_t lower_bits(uint64_t x)
14 {
15         return (x & ((ISAAC64_SZ-1) << 3)) >>3;
16 }
17
18 /* Extract next ISAAC64_SZ_LOG bits (starting at bit ISAAC64_SZ_LOG+2). */
19 static inline uint32_t upper_bits(uint32_t y)
20 {
21         return (y >> (ISAAC64_SZ_LOG+3)) & (ISAAC64_SZ-1);
22 }
23
24 static void isaac64_update(isaac64_ctx *_ctx){
25   uint64_t *m;
26   uint64_t *r;
27   uint64_t  a;
28   uint64_t  b;
29   uint64_t  x;
30   uint64_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<ISAAC64_SZ/2;i++){
37     x=m[i];
38     a=~(a^a<<21)+m[i+ISAAC64_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>>5)+m[i+ISAAC64_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<<12)+m[i+ISAAC64_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>>33)+m[i+ISAAC64_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=ISAAC64_SZ/2;i<ISAAC64_SZ;i++){
55     x=m[i];
56     a=~(a^a<<21)+m[i-ISAAC64_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>>5)+m[i-ISAAC64_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<<12)+m[i-ISAAC64_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>>33)+m[i-ISAAC64_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=ISAAC64_SZ;
75 }
76
77 static void isaac64_mix(uint64_t _x[8]){
78   static const unsigned char SHIFT[8]={9,9,23,15,14,20,17,14};
79   int i;
80   for(i=0;i<8;i++){
81     _x[i]-=_x[(i+4)&7];
82     _x[(i+5)&7]^=_x[(i+7)&7]>>SHIFT[i];
83     _x[(i+7)&7]+=_x[i];
84     i++;
85     _x[i]-=_x[(i+4)&7];
86     _x[(i+5)&7]^=_x[(i+7)&7]<<SHIFT[i];
87     _x[(i+7)&7]+=_x[i];
88   }
89 }
90
91
92 void isaac64_init(isaac64_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   isaac64_reseed(_ctx,_seed,_nseed);
96 }
97
98 void isaac64_reseed(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed){
99   uint64_t *m;
100   uint64_t *r;
101   uint64_t  x[8];
102   int       i;
103   int       j;
104   m=_ctx->m;
105   r=_ctx->r;
106   if(_nseed>ISAAC64_SEED_SZ_MAX)_nseed=ISAAC64_SEED_SZ_MAX;
107   for(i=0;i<_nseed>>3;i++){
108     r[i]^=(uint64_t)_seed[i<<3|7]<<56|(uint64_t)_seed[i<<3|6]<<48|
109      (uint64_t)_seed[i<<3|5]<<40|(uint64_t)_seed[i<<3|4]<<32|
110      (uint64_t)_seed[i<<3|3]<<24|(uint64_t)_seed[i<<3|2]<<16|
111      (uint64_t)_seed[i<<3|1]<<8|_seed[i<<3];
112   }
113   _nseed-=i<<3;
114   if(_nseed>0){
115     uint64_t ri;
116     ri=_seed[i<<3];
117     for(j=1;j<_nseed;j++)ri|=(uint64_t)_seed[i<<3|j]<<(j<<3);
118     r[i++]^=ri;
119   }
120   x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=(uint64_t)0x9E3779B97F4A7C13ULL;
121   for(i=0;i<4;i++)isaac64_mix(x);
122   for(i=0;i<ISAAC64_SZ;i+=8){
123     for(j=0;j<8;j++)x[j]+=r[i+j];
124     isaac64_mix(x);
125     memcpy(m+i,x,sizeof(x));
126   }
127   for(i=0;i<ISAAC64_SZ;i+=8){
128     for(j=0;j<8;j++)x[j]+=m[i+j];
129     isaac64_mix(x);
130     memcpy(m+i,x,sizeof(x));
131   }
132   isaac64_update(_ctx);
133 }
134
135 uint64_t isaac64_next_uint64(isaac64_ctx *_ctx){
136   if(!_ctx->n)isaac64_update(_ctx);
137   return _ctx->r[--_ctx->n];
138 }
139
140 uint64_t isaac64_next_uint(isaac64_ctx *_ctx,uint64_t _n){
141   uint64_t r;
142   uint64_t v;
143   uint64_t d;
144   do{
145     r=isaac64_next_uint64(_ctx);
146     v=r%_n;
147     d=r-v;
148   }
149   while(((d+_n-1)&ISAAC64_MASK)<d);
150   return v;
151 }
152
153 /*Returns a uniform random float.
154   The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
155   _bits: An initial set of random bits.
156   _base: This should be -(the number of bits in _bits), up to -64.
157   Return: A float uniformly distributed between 0 (inclusive) and 1
158            (exclusive).
159           The average value was measured over 2**32 samples to be
160            0.499991407275206357.*/
161 static float isaac64_float_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
162   float ret;
163   int   nbits_needed;
164   while(!_bits){
165     if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
166     _base-=64;
167     _bits=isaac64_next_uint64(_ctx);
168   }
169   nbits_needed=FLT_MANT_DIG-ILOGNZ_64(_bits);
170 #if FLT_MANT_DIG>64
171   ret=ldexpf((float)_bits,_base);
172 # if FLT_MANT_DIG>129
173   while(64-nbits_needed<0){
174 # else
175   if(64-nbits_needed<0){
176 # endif
177     _base-=64;
178     nbits_needed-=64;
179     ret+=ldexpf((float)isaac64_next_uint64(_ctx),_base);
180   }
181   _bits=isaac64_next_uint64(_ctx)>>(64-nbits_needed);
182   ret+=ldexpf((float)_bits,_base-nbits_needed);
183 #else
184   if(nbits_needed>0){
185     _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>(64-nbits_needed);
186   }
187 # if FLT_MANT_DIG<64
188   else _bits>>=-nbits_needed;
189 # endif
190   ret=ldexpf((float)_bits,_base-nbits_needed);
191 #endif
192   return ret;
193 }
194
195 float isaac64_next_float(isaac64_ctx *_ctx){
196   return isaac64_float_bits(_ctx,0,0);
197 }
198
199 float isaac64_next_signed_float(isaac64_ctx *_ctx){
200   uint64_t bits;
201   bits=isaac64_next_uint64(_ctx);
202   return (1|-((int)bits&1))*isaac64_float_bits(_ctx,bits>>1,-63);
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 -64.
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.499990992392019273.*/
212 static double isaac64_double_bits(isaac64_ctx *_ctx,uint64_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-=64;
218     _bits=isaac64_next_uint64(_ctx);
219   }
220   nbits_needed=DBL_MANT_DIG-ILOGNZ_64(_bits);
221 #if DBL_MANT_DIG>64
222   ret=ldexp((double)_bits,_base);
223 # if DBL_MANT_DIG>129
224   while(64-nbits_needed<0){
225 # else
226   if(64-nbits_needed<0){
227 # endif
228     _base-=64;
229     nbits_needed-=64;
230     ret+=ldexp((double)isaac64_next_uint64(_ctx),_base);
231   }
232   _bits=isaac64_next_uint64(_ctx)>>(64-nbits_needed);
233   ret+=ldexp((double)_bits,_base-nbits_needed);
234 #else
235   if(nbits_needed>0){
236     _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>(64-nbits_needed);
237   }
238 # if DBL_MANT_DIG<64
239   else _bits>>=-nbits_needed;
240 # endif
241   ret=ldexp((double)_bits,_base-nbits_needed);
242 #endif
243   return ret;
244 }
245
246 double isaac64_next_double(isaac64_ctx *_ctx){
247   return isaac64_double_bits(_ctx,0,0);
248 }
249
250 double isaac64_next_signed_double(isaac64_ctx *_ctx){
251   uint64_t bits;
252   bits=isaac64_next_uint64(_ctx);
253   return (1|-((int)bits&1))*isaac64_double_bits(_ctx,bits>>1,-63);
254 }