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