git.ozlabs.org Git - ccan/blob - 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 }