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 #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 ISAAC_MASK        (0xFFFFFFFFU)
  17
  18
  19
  20 static void isaac_update(isaac_ctx *_ctx){
  21   uint32_t *m;
  22   uint32_t *r;
  23   uint32_t  a;
  24   uint32_t  b;
  25   uint32_t  x;
  26   uint32_t  y;
  27   int       i;
  28   m=_ctx->m;
  29   r=_ctx->r;
  30   a=_ctx->a;
  31   b=_ctx->b+(++_ctx->c)&ISAAC_MASK;
  32   for(i=0;i<ISAAC_SZ/2;i++){
  33     x=m[i];
  34     a=(a^a<<13)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
  35     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  36     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  37     x=m[++i];
  38     a=(a^a>>6)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
  39     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  40     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  41     x=m[++i];
  42     a=(a^a<<2)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
  43     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  44     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  45     x=m[++i];
  46     a=(a^a>>16)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
  47     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  48     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  49   }
  50   for(i=ISAAC_SZ/2;i<ISAAC_SZ;i++){
  51     x=m[i];
  52     a=(a^a<<13)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
  53     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  54     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  55     x=m[++i];
  56     a=(a^a>>6)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
  57     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  58     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  59     x=m[++i];
  60     a=(a^a<<2)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
  61     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  62     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  63     x=m[++i];
  64     a=(a^a>>16)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
  65     m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
  66     r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
  67   }
  68   _ctx->b=b;
  69   _ctx->a=a;
  70   _ctx->n=ISAAC_SZ;
  71 }
  72
  73 static void isaac_mix(uint32_t _x[8]){
  74   static const unsigned char SHIFT[8]={11,2,8,16,10,4,8,9};
  75   int i;
  76   for(i=0;i<8;i++){
  77     _x[i]^=_x[i+1&7]<<SHIFT[i];
  78     _x[i+3&7]+=_x[i];
  79     _x[i+1&7]+=_x[i+2&7];
  80     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   }
  85 }
  86
  87
  88 void isaac_init(isaac_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   isaac_reseed(_ctx,_seed,_nseed);
  92 }
  93
  94 void isaac_reseed(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
  95   uint32_t *m;
  96   uint32_t *r;
  97   uint32_t  x[8];
  98   int       i;
  99   int       j;
 100   m=_ctx->m;
 101   r=_ctx->r;
 102   if(_nseed>ISAAC_SEED_SZ_MAX)_nseed=ISAAC_SEED_SZ_MAX;
 103   for(i=0;i<_nseed>>2;i++){
 104     r[i]^=(uint32_t)_seed[i<<2|3]<<24|(uint32_t)_seed[i<<2|2]<<16|
 105      (uint32_t)_seed[i<<2|1]<<8|_seed[i<<2];
 106   }
 107   _nseed-=i<<2;
 108   if(_nseed>0){
 109     uint32_t ri;
 110     ri=_seed[i<<2];
 111     for(j=1;j<_nseed;j++)ri|=(uint32_t)_seed[i<<2|j]<<(j<<3);
 112     r[i++]^=ri;
 113   }
 114   x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=0x9E3779B9U;
 115   for(i=0;i<4;i++)isaac_mix(x);
 116   for(i=0;i<ISAAC_SZ;i+=8){
 117     for(j=0;j<8;j++)x[j]+=r[i+j];
 118     isaac_mix(x);
 119     memcpy(m+i,x,sizeof(x));
 120   }
 121   for(i=0;i<ISAAC_SZ;i+=8){
 122     for(j=0;j<8;j++)x[j]+=m[i+j];
 123     isaac_mix(x);
 124     memcpy(m+i,x,sizeof(x));
 125   }
 126   isaac_update(_ctx);
 127 }
 128
 129 uint32_t isaac_next_uint32(isaac_ctx *_ctx){
 130   if(!_ctx->n)isaac_update(_ctx);
 131   return _ctx->r[--_ctx->n];
 132 }
 133
 134 uint32_t isaac_next_uint(isaac_ctx *_ctx,uint32_t _n){
 135   uint32_t r;
 136   uint32_t v;
 137   uint32_t d;
 138   do{
 139     r=isaac_next_uint32(_ctx);
 140     v=r%_n;
 141     d=r-v;
 142   }
 143   while((d+_n-1&ISAAC_MASK)<d);
 144   return v;
 145 }
 146
 147 /*Returns a uniform random float.
 148   The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
 149   _bits: An initial set of random bits.
 150   _base: This should be -(the number of bits in _bits), up to -32.
 151   Return: A float uniformly distributed between 0 (inclusive) and 1
 152            (exclusive).
 153           The average value was measured over 2**32 samples to be
 154            0.50000037448772916.*/
 155 static float isaac_float_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
 156   float ret;
 157   int   nbits_needed;
 158   while(!_bits){
 159     if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
 160     _base-=32;
 161     _bits=isaac_next_uint32(_ctx);
 162   }
 163   /*Note: This could also be determined with frexp(), for a slightly more
 164      portable solution, but that takes twice as long, and one has to worry
 165      about rounding effects, which can over-estimate the exponent when given
 166      FLT_MANT_DIG+1 consecutive one bits.
 167     Even the fallback C implementation of ILOGNZ_32() yields an implementation
 168      25% faster than the frexp() method.*/
 169   nbits_needed=FLT_MANT_DIG-ILOGNZ_32(_bits);
 170 #if FLT_MANT_DIG>32
 171   ret=ldexpf((float)_bits,_base);
 172 # if FLT_MANT_DIG>65
 173   while(32-nbits_needed<0){
 174 # else
 175   if(32-nbits_needed<0){
 176 # endif
 177     _base-=32;
 178     nbits_needed-=32;
 179     ret+=ldexpf((float)isaac_next_uint32(_ctx),_base);
 180   }
 181   _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
 182   ret+=ldexpf((float)_bits,_base-nbits_needed);
 183 #else
 184   if(nbits_needed>0){
 185     _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
 186   }
 187 # if FLT_MANT_DIG<32
 188   else _bits>>=-nbits_needed;
 189 # endif
 190   ret=ldexpf((float)_bits,_base-nbits_needed);
 191 #endif
 192   return ret;
 193 }
 194
 195 float isaac_next_float(isaac_ctx *_ctx){
 196   return isaac_float_bits(_ctx,0,0);
 197 }
 198
 199 float isaac_next_signed_float(isaac_ctx *_ctx){
 200   uint32_t bits;
 201   bits=isaac_next_uint32(_ctx);
 202   return (1|-((int)bits&1))*isaac_float_bits(_ctx,bits>>1,-31);
 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 -32.
 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.500006289408060911*/
 212 static double isaac_double_bits(isaac_ctx *_ctx,uint32_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-=32;
 218     _bits=isaac_next_uint32(_ctx);
 219   }
 220   nbits_needed=DBL_MANT_DIG-ILOGNZ_32(_bits);
 221 #if DBL_MANT_DIG>32
 222   ret=ldexp((double)_bits,_base);
 223 # if DBL_MANT_DIG>65
 224   while(32-nbits_needed<0){
 225 # else
 226   if(32-nbits_needed<0){
 227 # endif
 228     _base-=32;
 229     nbits_needed-=32;
 230     ret+=ldexp((double)isaac_next_uint32(_ctx),_base);
 231   }
 232   _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
 233   ret+=ldexp((double)_bits,_base-nbits_needed);
 234 #else
 235   if(nbits_needed>0){
 236     _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
 237   }
 238 # if DBL_MANT_DIG<32
 239   else _bits>>=-nbits_needed;
 240 # endif
 241   ret=ldexp((double)_bits,exp-DBL_MANT_DIG);
 242 #endif
 243   return ret;
 244 }
 245
 246 double isaac_next_double(isaac_ctx *_ctx){
 247   return isaac_double_bits(_ctx,0,0);
 248 }
 249
 250 double isaac_next_signed_double(isaac_ctx *_ctx){
 251   uint32_t bits;
 252   bits=isaac_next_uint32(_ctx);
 253   return (1|-((int)bits&1))*isaac_double_bits(_ctx,bits>>1,-31);
 254 }