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