git.ozlabs.org Git - ccan/blob - ccan/isaac/isaac.c

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