Tim's ISAAC module.
[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 #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 ISAAC64_MASK ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
17
18
19
20 static void isaac64_update(isaac64_ctx *_ctx){
21   uint64_t *m;
22   uint64_t *r;
23   uint64_t  a;
24   uint64_t  b;
25   uint64_t  x;
26   uint64_t  y;
27   int       i;
28   m=_ctx->m;
29   r=_ctx->r;
30   a=_ctx->a;
31   b=_ctx->b+(++_ctx->c)&ISAAC64_MASK;
32   for(i=0;i<ISAAC64_SZ/2;i++){
33     x=m[i];
34     a=~(a^a<<21)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
35     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
36     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
37     x=m[++i];
38     a=(a^a>>5)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
39     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
40     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
41     x=m[++i];
42     a=(a^a<<12)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
43     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
44     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
45     x=m[++i];
46     a=(a^a>>33)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
47     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
48     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
49   }
50   for(i=ISAAC64_SZ/2;i<ISAAC64_SZ;i++){
51     x=m[i];
52     a=~(a^a<<21)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
53     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
54     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
55     x=m[++i];
56     a=(a^a>>5)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
57     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
58     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
59     x=m[++i];
60     a=(a^a<<12)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
61     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
62     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
63     x=m[++i];
64     a=(a^a>>33)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
65     m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
66     r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
67   }
68   _ctx->b=b;
69   _ctx->a=a;
70   _ctx->n=ISAAC64_SZ;
71 }
72
73 static void isaac64_mix(uint64_t _x[8]){
74   static const unsigned char SHIFT[8]={9,9,23,15,14,20,17,14};
75   int i;
76   for(i=0;i<8;i++){
77     _x[i]-=_x[i+4&7];
78     _x[i+5&7]^=_x[i+7&7]>>SHIFT[i];
79     _x[i+7&7]+=_x[i];
80     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   }
85 }
86
87
88 void isaac64_init(isaac64_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   isaac64_reseed(_ctx,_seed,_nseed);
92 }
93
94 void isaac64_reseed(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed){
95   uint64_t *m;
96   uint64_t *r;
97   uint64_t  x[8];
98   int       i;
99   int       j;
100   m=_ctx->m;
101   r=_ctx->r;
102   if(_nseed>ISAAC64_SEED_SZ_MAX)_nseed=ISAAC64_SEED_SZ_MAX;
103   for(i=0;i<_nseed>>3;i++){
104     r[i]^=(uint64_t)_seed[i<<3|7]<<56|(uint64_t)_seed[i<<3|6]<<48|
105      (uint64_t)_seed[i<<3|5]<<40|(uint64_t)_seed[i<<3|4]<<32|
106      (uint64_t)_seed[i<<3|3]<<24|(uint64_t)_seed[i<<3|2]<<16|
107      (uint64_t)_seed[i<<3|1]<<8|_seed[i<<3];
108   }
109   _nseed-=i<<3;
110   if(_nseed>0){
111     uint64_t ri;
112     ri=_seed[i<<3];
113     for(j=1;j<_nseed;j++)ri|=(uint64_t)_seed[i<<3|j]<<(j<<3);
114     r[i++]^=ri;
115   }
116   x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=(uint64_t)0x9E3779B97F4A7C13ULL;
117   for(i=0;i<4;i++)isaac64_mix(x);
118   for(i=0;i<ISAAC64_SZ;i+=8){
119     for(j=0;j<8;j++)x[j]+=r[i+j];
120     isaac64_mix(x);
121     memcpy(m+i,x,sizeof(x));
122   }
123   for(i=0;i<ISAAC64_SZ;i+=8){
124     for(j=0;j<8;j++)x[j]+=m[i+j];
125     isaac64_mix(x);
126     memcpy(m+i,x,sizeof(x));
127   }
128   isaac64_update(_ctx);
129 }
130
131 uint64_t isaac64_next_uint64(isaac64_ctx *_ctx){
132   if(!_ctx->n)isaac64_update(_ctx);
133   return _ctx->r[--_ctx->n];
134 }
135
136 uint64_t isaac64_next_uint(isaac64_ctx *_ctx,uint64_t _n){
137   uint64_t r;
138   uint64_t v;
139   uint64_t d;
140   do{
141     r=isaac64_next_uint64(_ctx);
142     v=r%_n;
143     d=r-v;
144   }
145   while((d+_n-1&ISAAC64_MASK)<d);
146   return v;
147 }
148
149 /*Returns a uniform random float.
150   The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
151   _bits: An initial set of random bits.
152   _base: This should be -(the number of bits in _bits), up to -64.
153   Return: A float uniformly distributed between 0 (inclusive) and 1
154            (exclusive).
155           The average value was measured over 2**32 samples to be
156            0.499991407275206357.*/
157 static float isaac64_float_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
158   float ret;
159   int   nbits_needed;
160   while(!_bits){
161     if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
162     _base-=64;
163     _bits=isaac64_next_uint64(_ctx);
164   }
165   nbits_needed=FLT_MANT_DIG-ILOGNZ_64(_bits);
166 #if FLT_MANT_DIG>64
167   ret=ldexpf((float)_bits,_base);
168 # if FLT_MANT_DIG>129
169   while(64-nbits_needed<0){
170 # else
171   if(64-nbits_needed<0){
172 # endif
173     _base-=64;
174     nbits_needed-=64;
175     ret+=ldexpf((float)isaac64_next_uint64(_ctx),_base);
176   }
177   _bits=isaac64_next_uint64(_ctx)>>64-nbits_needed;
178   ret+=ldexpf((float)_bits,_base-nbits_needed);
179 #else
180   if(nbits_needed>0){
181     _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>64-nbits_needed;
182   }
183 # if FLT_MANT_DIG<64
184   else _bits>>=-nbits_needed;
185 # endif
186   ret=ldexpf((float)_bits,_base-nbits_needed);
187 #endif
188   return ret;
189 }
190
191 float isaac64_next_float(isaac64_ctx *_ctx){
192   return isaac64_float_bits(_ctx,0,0);
193 }
194
195 float isaac64_next_signed_float(isaac64_ctx *_ctx){
196   uint64_t bits;
197   bits=isaac64_next_uint64(_ctx);
198   return (1|-((int)bits&1))*isaac64_float_bits(_ctx,bits>>1,-63);
199 }
200
201 /*Returns a uniform random double.
202   _bits: An initial set of random bits.
203   _base: This should be -(the number of bits in _bits), up to -64.
204   Return: A double uniformly distributed between 0 (inclusive) and 1
205            (exclusive).
206           The average value was measured over 2**32 samples to be
207            0.499990992392019273.*/
208 static double isaac64_double_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
209   double ret;
210   int    nbits_needed;
211   while(!_bits){
212     if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
213     _base-=64;
214     _bits=isaac64_next_uint64(_ctx);
215   }
216   nbits_needed=DBL_MANT_DIG-ILOGNZ_64(_bits);
217 #if DBL_MANT_DIG>64
218   ret=ldexp((double)_bits,_base);
219 # if DBL_MANT_DIG>129
220   while(64-nbits_needed<0){
221 # else
222   if(64-nbits_needed<0){
223 # endif
224     _base-=64;
225     nbits_needed-=64;
226     ret+=ldexp((double)isaac64_next_uint64(_ctx),_base);
227   }
228   _bits=isaac64_next_uint64(_ctx)>>64-nbits_needed;
229   ret+=ldexp((double)_bits,_base-nbits_needed);
230 #else
231   if(nbits_needed>0){
232     _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>64-nbits_needed;
233   }
234 # if DBL_MANT_DIG<64
235   else _bits>>=-nbits_needed;
236 # endif
237   ret=ldexp((double)_bits,_base-nbits_needed);
238 #endif
239   return ret;
240 }
241
242 double isaac64_next_double(isaac64_ctx *_ctx){
243   return isaac64_double_bits(_ctx,0,0);
244 }
245
246 double isaac64_next_signed_double(isaac64_ctx *_ctx){
247   uint64_t bits;
248   bits=isaac64_next_uint64(_ctx);
249   return (1|-((int)bits&1))*isaac64_double_bits(_ctx,bits>>1,-63);
250 }