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