Tim's ISAAC module.

ccan/isaac/_info.c

+/**
+ * isaac - A fast, high-quality pseudo-random number generator.
+ *
+ * ISAAC (Indirect, Shift, Accumulate, Add, and Count) is the most advanced of
+ *  a series of pseudo-random number generators designed by Robert J. Jenkins
+ *  Jr. in 1996: http://www.burtleburtle.net/bob/rand/isaac.html
+ * To quote:
+ *   No efficient method is known for deducing their internal states.
+ *   ISAAC requires an amortized 18.75 instructions to produce a 32-bit value.
+ *   There are no cycles in ISAAC shorter than 2**40 values.
+ *   The expected cycle length is 2**8295 values.
+ *   ...
+ *   ISAAC-64 generates a different sequence than ISAAC, but it uses the same
+ *    principles.
+ *   It uses 64-bit arithmetic.
+ *   It generates a 64-bit result every 19 instructions.
+ *   All cycles are at least 2**72 values, and the average cycle length is
+ *    2**16583.
+ * An additional, important comment from Bob Jenkins in 2006:
+ *   Seeding a random number generator is essentially the same problem as
+ *    encrypting the seed with a block cipher.
+ *   ISAAC should be initialized with the encryption of the seed by some
+ *    secure cipher.
+ *   I've provided a seeding routine in my implementations, which nobody has
+ *    broken so far, but I have less faith in that initialization routine than
+ *    I have in ISAAC.
+ *
+ * A number of attacks on ISAAC have been published.
+ * [Pudo01] can recover the entire internal state and has expected running time
+ *  less than the square root of the number of states, or 2**4121 (4.67E+1240).
+ * [Auma06] reveals a large set of weak states, consisting of those for which
+ *  the first value is repeated one or more times elsewhere in the state
+ *  vector.
+ * These induce a bias in the output relative to the repeated value.
+ * The seed values used as input below are scrambled before being used, so any
+ *  duplicates in them do not imply duplicates in the resulting internal state,
+ *  however the chances of some duplicate existing elsewhere in a random state
+ *  are just over 255/2**32, or merely 1 in 16 million.
+ * Such states are, of course, much rarer in ISAAC-64.
+ * It is not clear if an attacker can tell from just the output if ISAAC is in
+ *  a weak state, or deduce the full internal state in any case except that
+ *  where all or almost all of the entries in the state vector are identical.
+ *   @MISC{Pudo01,
+ *     author="Marina Pudovkina",
+ *     title="A Known Plaintext Attack on the {ISAAC} Keystream Generator",
+ *     howpublished="Cryptology ePrint Archive, Report 2001/049",
+ *     year=2001,
+ *     note="\url{http://eprint.iacr.org/2001/049}",
+ *   }
+ *   @MISC{Auma06,
+ *     author="Jean-Philippe Aumasson",
+ *     title="On the Pseudo-Random Generator {ISAAC}",
+ *     howpublished="Cryptology ePrint Archive, Report 2006/438",
+ *     year=2006,
+ *     note="\url{http://eprint.iacr.org/2006/438}",
+ *   }
+ *
+ * Even if one does not trust the security of this PRNG (and, without a good
+ *  source of entropy to seed it, one should not), ISAAC is an excellent source
+ *  of high-quality random numbers for Monte Carlo simulations, etc.
+ * It is the fastest 32-bit generator among all of those that pass the
+ *  statistical tests in the recent survey
+ *  http://www.iro.umontreal.ca/~simardr/testu01/tu01.html, with the exception
+ *  of Marsa-LFIB4, and it is quite competitive on 64-bit archtectures.
+ * Unlike Marsa-LFIB4 (and all other LFib generators), there are no linear
+ *  dependencies between successive values, and unlike many generators found in
+ *  libc implementations, there are no small periods in the least significant
+ *  bits, or seeds which lead to very small periods in general.
+ *
+ * Example:
+ *  #include <stdio.h>
+ *  #include <time.h>
+ *  #include <ccan/isaac/isaac.h>
+ *
+ *  int main(void){
+ *    static const char *CHEESE[3]={"Cheddar","Provolone","Camembert"};
+ *    isaac_ctx     isaac;
+ *    unsigned char seed[8];
+ *    time_t        now;
+ *    int           i;
+ *    //N.B.: time() is not a good source of entropy.
+ *    //Do not use it for cryptogrpahic purposes.
+ *    time(&now);
+ *    //Print it out so we can reproduce problems if needed.
+ *    printf("Seed: 0x%016llX\n",(long long)now);
+ *    //And convert the time to a byte array so that we can reproduce the same
+ *    // seed on platforms with different endianesses.
+ *    for(i=0;i<8;i++){
+ *      seed[i]=(unsigned char)(now&0xFF);
+ *      now>>=8;
+ *    }
+ *    isaac_init(&isaac,seed,8);
+ *    printf("0x%08lX\n",(long)isaac_next_uint32(&isaac));
+ *    printf("%s\n",CHEESE[isaac_next_uint(&isaac,3)]);
+ *    printf("%0.8G\n",isaac_next_float(&isaac));
+ *    printf("%0.8G\n",isaac_next_signed_float(&isaac));
+ *    printf("%0.18G\n",isaac_next_double(&isaac));
+ *    printf("%0.18G\n",isaac_next_signed_double(&isaac));
+ *    return 0;
+ *  }
+ *
+ * License: Public Domain
+ */
+#include <string.h>
+#include "config.h"
+
+int main(int _argc,const char *_argv[]){
+  /*Expect exactly one argument.*/
+  if(_argc!=2)return 1;
+  if(strcmp(_argv[1],"depends")==0){
+    /*PRINTF-CCAN-PACKAGES-YOU-NEED-ONE-PER-LINE-IF-ANY*/
+    printf("ccan/ilog\n");
+    return 0;
+  }
+  return 1;
+}

ccan/isaac/isaac.c

+/*Written by Timothy B. Terriberry (tterribe@xiph.org) 1999-2009 public domain.
+  Based on the public domain implementation by Robert J. Jenkins Jr.*/
+#include <float.h>
+#include <math.h>
+#include <string.h>
+#include <ccan/ilog/ilog.h>
+#include "isaac.h"
+#if defined(__GNUC_PREREQ)
+# if __GNUC_PREREQ(4,2)
+#  pragma GCC diagnostic ignored "-Wparentheses"
+# endif
+#endif
+
+
+
+#define ISAAC_MASK        (0xFFFFFFFFU)
+
+
+
+static void isaac_update(isaac_ctx *_ctx){
+  uint32_t *m;
+  uint32_t *r;
+  uint32_t  a;
+  uint32_t  b;
+  uint32_t  x;
+  uint32_t  y;
+  int       i;
+  m=_ctx->m;
+  r=_ctx->r;
+  a=_ctx->a;
+  b=_ctx->b+(++_ctx->c)&ISAAC_MASK;
+  for(i=0;i<ISAAC_SZ/2;i++){
+    x=m[i];
+    a=(a^a<<13)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+    x=m[++i];
+    a=(a^a>>6)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+    x=m[++i];
+    a=(a^a<<2)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+    x=m[++i];
+    a=(a^a>>16)+m[i+ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+  }
+  for(i=ISAAC_SZ/2;i<ISAAC_SZ;i++){
+    x=m[i];
+    a=(a^a<<13)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+    x=m[++i];
+    a=(a^a>>6)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+    x=m[++i];
+    a=(a^a<<2)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+    x=m[++i];
+    a=(a^a>>16)+m[i-ISAAC_SZ/2]&ISAAC_MASK;
+    m[i]=y=m[(x&ISAAC_SZ-1<<2)>>2]+a+b&ISAAC_MASK;
+    r[i]=b=m[y>>ISAAC_SZ_LOG+2&ISAAC_SZ-1]+x&ISAAC_MASK;
+  }
+  _ctx->b=b;
+  _ctx->a=a;
+  _ctx->n=ISAAC_SZ;
+}
+
+static void isaac_mix(uint32_t _x[8]){
+  static const unsigned char SHIFT[8]={11,2,8,16,10,4,8,9};
+  int i;
+  for(i=0;i<8;i++){
+    _x[i]^=_x[i+1&7]<<SHIFT[i];
+    _x[i+3&7]+=_x[i];
+    _x[i+1&7]+=_x[i+2&7];
+    i++;
+    _x[i]^=_x[i+1&7]>>SHIFT[i];
+    _x[i+3&7]+=_x[i];
+    _x[i+1&7]+=_x[i+2&7];
+  }
+}
+
+
+void isaac_init(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
+  _ctx->a=_ctx->b=_ctx->c=0;
+  memset(_ctx->r,0,sizeof(_ctx->r));
+  isaac_reseed(_ctx,_seed,_nseed);
+}
+
+void isaac_reseed(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed){
+  uint32_t *m;
+  uint32_t *r;
+  uint32_t  x[8];
+  int       i;
+  int       j;
+  m=_ctx->m;
+  r=_ctx->r;
+  if(_nseed>ISAAC_SEED_SZ_MAX)_nseed=ISAAC_SEED_SZ_MAX;
+  for(i=0;i<_nseed>>2;i++){
+    r[i]^=(uint32_t)_seed[i<<2|3]<<24|(uint32_t)_seed[i<<2|2]<<16|
+     (uint32_t)_seed[i<<2|1]<<8|_seed[i<<2];
+  }
+  _nseed-=i<<2;
+  if(_nseed>0){
+    uint32_t ri;
+    ri=_seed[i<<2];
+    for(j=1;j<_nseed;j++)ri|=(uint32_t)_seed[i<<2|j]<<(j<<3);
+    r[i++]^=ri;
+  }
+  x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=0x9E3779B9U;
+  for(i=0;i<4;i++)isaac_mix(x);
+  for(i=0;i<ISAAC_SZ;i+=8){
+    for(j=0;j<8;j++)x[j]+=r[i+j];
+    isaac_mix(x);
+    memcpy(m+i,x,sizeof(x));
+  }
+  for(i=0;i<ISAAC_SZ;i+=8){
+    for(j=0;j<8;j++)x[j]+=m[i+j];
+    isaac_mix(x);
+    memcpy(m+i,x,sizeof(x));
+  }
+  isaac_update(_ctx);
+}
+
+uint32_t isaac_next_uint32(isaac_ctx *_ctx){
+  if(!_ctx->n)isaac_update(_ctx);
+  return _ctx->r[--_ctx->n];
+}
+
+uint32_t isaac_next_uint(isaac_ctx *_ctx,uint32_t _n){
+  uint32_t r;
+  uint32_t v;
+  uint32_t d;
+  do{
+    r=isaac_next_uint32(_ctx);
+    v=r%_n;
+    d=r-v;
+  }
+  while((d+_n-1&ISAAC_MASK)<d);
+  return v;
+}
+
+/*Returns a uniform random float.
+  The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
+  _bits: An initial set of random bits.
+  _base: This should be -(the number of bits in _bits), up to -32.
+  Return: A float uniformly distributed between 0 (inclusive) and 1
+           (exclusive).
+          The average value was measured over 2**32 samples to be
+           0.50000037448772916.*/
+static float isaac_float_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
+  float ret;
+  int   nbits_needed;
+  while(!_bits){
+    if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
+    _base-=32;
+    _bits=isaac_next_uint32(_ctx);
+  }
+  /*Note: This could also be determined with frexp(), for a slightly more
+     portable solution, but that takes twice as long, and one has to worry
+     about rounding effects, which can over-estimate the exponent when given
+     FLT_MANT_DIG+1 consecutive one bits.
+    Even the fallback C implementation of ILOGNZ_32() yields an implementation
+     25% faster than the frexp() method.*/
+  nbits_needed=FLT_MANT_DIG-ILOGNZ_32(_bits);
+#if FLT_MANT_DIG>32
+  ret=ldexpf((float)_bits,_base);
+# if FLT_MANT_DIG>65
+  while(32-nbits_needed<0){
+# else
+  if(32-nbits_needed<0){
+# endif
+    _base-=32;
+    nbits_needed-=32;
+    ret+=ldexpf((float)isaac_next_uint32(_ctx),_base);
+  }
+  _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
+  ret+=ldexpf((float)_bits,_base-nbits_needed);
+#else
+  if(nbits_needed>0){
+    _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
+  }
+# if FLT_MANT_DIG<32
+  else _bits>>=-nbits_needed;
+# endif
+  ret=ldexpf((float)_bits,_base-nbits_needed);
+#endif
+  return ret;
+}
+
+float isaac_next_float(isaac_ctx *_ctx){
+  return isaac_float_bits(_ctx,0,0);
+}
+
+float isaac_next_signed_float(isaac_ctx *_ctx){
+  uint32_t bits;
+  bits=isaac_next_uint32(_ctx);
+  return (1|-((int)bits&1))*isaac_float_bits(_ctx,bits>>1,-31);
+}
+
+/*Returns a uniform random double.
+  _bits: An initial set of random bits.
+  _base: This should be -(the number of bits in _bits), up to -32.
+  Return: A double uniformly distributed between 0 (inclusive) and 1
+           (exclusive).
+          The average value was measured over 2**32 samples to be
+           0.500006289408060911*/
+static double isaac_double_bits(isaac_ctx *_ctx,uint32_t _bits,int _base){
+  double ret;
+  int    nbits_needed;
+  while(!_bits){
+    if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
+    _base-=32;
+    _bits=isaac_next_uint32(_ctx);
+  }
+  nbits_needed=DBL_MANT_DIG-ILOGNZ_32(_bits);
+#if DBL_MANT_DIG>32
+  ret=ldexp((double)_bits,_base);
+# if DBL_MANT_DIG>65
+  while(32-nbits_needed<0){
+# else
+  if(32-nbits_needed<0){
+# endif
+    _base-=32;
+    nbits_needed-=32;
+    ret+=ldexp((double)isaac_next_uint32(_ctx),_base);
+  }
+  _bits=isaac_next_uint32(_ctx)>>32-nbits_needed;
+  ret+=ldexp((double)_bits,_base-nbits_needed);
+#else
+  if(nbits_needed>0){
+    _bits=_bits<<nbits_needed|isaac_next_uint32(_ctx)>>32-nbits_needed;
+  }
+# if DBL_MANT_DIG<32
+  else _bits>>=-nbits_needed;
+# endif
+  ret=ldexp((double)_bits,exp-DBL_MANT_DIG);
+#endif
+  return ret;
+}
+
+double isaac_next_double(isaac_ctx *_ctx){
+  return isaac_double_bits(_ctx,0,0);
+}
+
+double isaac_next_signed_double(isaac_ctx *_ctx){
+  uint32_t bits;
+  bits=isaac_next_uint32(_ctx);
+  return (1|-((int)bits&1))*isaac_double_bits(_ctx,bits>>1,-31);
+}

ccan/isaac/isaac.h

+#if !defined(_isaac_H)
+# define _isaac_H (1)
+# include <stdint.h>
+
+
+
+typedef struct isaac_ctx isaac_ctx;
+
+
+
+/*This value may be lowered to reduce memory usage on embedded platforms, at
+   the cost of reducing security and increasing bias.
+  Quoting Bob Jenkins: "The current best guess is that bias is detectable after
+   2**37 values for [ISAAC_SZ_LOG]=3, 2**45 for 4, 2**53 for 5, 2**61 for 6,
+   2**69 for 7, and 2**77 values for [ISAAC_SZ_LOG]=8."*/
+#define ISAAC_SZ_LOG      (8)
+#define ISAAC_SZ          (1<<ISAAC_SZ_LOG)
+#define ISAAC_SEED_SZ_MAX (ISAAC_SZ<<2)
+
+
+
+/*ISAAC is the most advanced of a series of pseudo-random number generators
+   designed by Robert J. Jenkins Jr. in 1996.
+  http://www.burtleburtle.net/bob/rand/isaac.html
+  To quote:
+    No efficient method is known for deducing their internal states.
+    ISAAC requires an amortized 18.75 instructions to produce a 32-bit value.
+    There are no cycles in ISAAC shorter than 2**40 values.
+    The expected cycle length is 2**8295 values.*/
+struct isaac_ctx{
+  unsigned n;
+  uint32_t r[ISAAC_SZ];
+  uint32_t m[ISAAC_SZ];
+  uint32_t a;
+  uint32_t b;
+  uint32_t c;
+};
+
+
+/**
+ * isaac_init - Initialize an instance of the ISAAC random number generator.
+ * @_ctx:   The instance to initialize.
+ * @_seed:  The specified seed bytes.
+ *          This may be NULL if _nseed is less than or equal to zero.
+ * @_nseed: The number of bytes to use for the seed.
+ *          If this is greater than ISAAC_SEED_SZ_MAX, the extra bytes are
+ *           ignored.
+ */
+void isaac_init(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed);
+
+/**
+ * isaac_reseed - Mix a new batch of entropy into the current state.
+ * To reset ISAAC to a known state, call isaac_init() again instead.
+ * @_ctx:   The instance to reseed.
+ * @_seed:  The specified seed bytes.
+ *          This may be NULL if _nseed is zero.
+ * @_nseed: The number of bytes to use for the seed.
+ *          If this is greater than ISAAC_SEED_SZ_MAX, the extra bytes are
+ *           ignored.
+ */
+void isaac_reseed(isaac_ctx *_ctx,const unsigned char *_seed,int _nseed);
+/**
+ * isaac_next_uint32 - Return the next random 32-bit value.
+ * @_ctx: The ISAAC instance to generate the value with.
+ */
+uint32_t isaac_next_uint32(isaac_ctx *_ctx);
+/**
+ * isaac_next_uint - Uniform random integer less than the given value.
+ * @_ctx: The ISAAC instance to generate the value with.
+ * @_n:   The upper bound on the range of numbers returned (not inclusive).
+ *        This must be greater than zero and less than 2**32.
+ *        To return integers in the full range 0...2**32-1, use
+ *         isaac_next_uint32() instead.
+ * Return: An integer uniformly distributed between 0 and _n-1 (inclusive).
+ */
+uint32_t isaac_next_uint(isaac_ctx *_ctx,uint32_t _n);
+/**
+ * isaac_next_float - Uniform random float in the range [0,1).
+ * @_ctx: The ISAAC instance to generate the value with.
+ * Returns a high-quality float uniformly distributed between 0 (inclusive)
+ *  and 1 (exclusive).
+ * All of the float's mantissa bits are random, e.g., the least significant bit
+ *  may still be non-zero even if the value is less than 0.5, and any
+ *  representable float in the range [0,1) has a chance to be returned, though
+ *  values very close to zero become increasingly unlikely.
+ * To generate cheaper float values that do not have these properties, use
+ *   ldexpf((float)isaac_next_uint32(_ctx),-32);
+ */
+float isaac_next_float(isaac_ctx *_ctx);
+/**
+ * isaac_next_signed_float - Uniform random float in the range (-1,1).
+ * @_ctx: The ISAAC instance to generate the value with.
+ * Returns a high-quality float uniformly distributed between -1 and 1
+ *  (exclusive).
+ * All of the float's mantissa bits are random, e.g., the least significant bit
+ *  may still be non-zero even if the magnitude is less than 0.5, and any
+ *  representable float in the range (-1,1) has a chance to be returned, though
+ *  values very close to zero become increasingly unlikely.
+ * To generate cheaper float values that do not have these properties, use
+ *   ldexpf((float)isaac_next_uint32(_ctx),-31)-1;
+ *  though this returns values in the range [-1,1).
+ */
+float isaac_next_signed_float(isaac_ctx *_ctx);
+/**
+ * isaac_next_double - Uniform random double in the range [0,1).
+ * @_ctx: The ISAAC instance to generate the value with.
+ * Returns a high-quality double uniformly distributed between 0 (inclusive)
+ *  and 1 (exclusive).
+ * All of the double's mantissa bits are random, e.g., the least significant
+ *  bit may still be non-zero even if the value is less than 0.5, and any
+ *  representable double in the range [0,1) has a chance to be returned, though
+ *  values very close to zero become increasingly unlikely.
+ * To generate cheaper double values that do not have these properties, use
+ *   ldexp((double)isaac_next_uint32(_ctx),-32);
+ */
+double isaac_next_double(isaac_ctx *_ctx);
+/**
+ * isaac_next_signed_double - Uniform random double in the range (-1,1).
+ * @_ctx: The ISAAC instance to generate the value with.
+ * Returns a high-quality double uniformly distributed between -1 and 1
+ *  (exclusive).
+ * All of the double's mantissa bits are random, e.g., the least significant
+ *  bit may still be non-zero even if the value is less than 0.5, and any
+ *  representable double in the range (-1,1) has a chance to be returned,
+ *  though values very close to zero become increasingly unlikely.
+ * To generate cheaper double values that do not have these properties, use
+ *   ldexp((double)isaac_next_uint32(_ctx),-31)-1;
+ *  though this returns values in the range [-1,1).
+ */
+double isaac_next_signed_double(isaac_ctx *_ctx);
+
+#endif

ccan/isaac/isaac64.c

+/*Written by Timothy B. Terriberry (tterribe@xiph.org) 1999-2009 public domain.
+  Based on the public domain ISAAC implementation by Robert J. Jenkins Jr.*/
+#include <float.h>
+#include <math.h>
+#include <string.h>
+#include <ccan/ilog/ilog.h>
+#include "isaac64.h"
+#if defined(__GNUC_PREREQ)
+# if __GNUC_PREREQ(4,2)
+#  pragma GCC diagnostic ignored "-Wparentheses"
+# endif
+#endif
+
+
+
+#define ISAAC64_MASK ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
+
+
+
+static void isaac64_update(isaac64_ctx *_ctx){
+  uint64_t *m;
+  uint64_t *r;
+  uint64_t  a;
+  uint64_t  b;
+  uint64_t  x;
+  uint64_t  y;
+  int       i;
+  m=_ctx->m;
+  r=_ctx->r;
+  a=_ctx->a;
+  b=_ctx->b+(++_ctx->c)&ISAAC64_MASK;
+  for(i=0;i<ISAAC64_SZ/2;i++){
+    x=m[i];
+    a=~(a^a<<21)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+    x=m[++i];
+    a=(a^a>>5)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+    x=m[++i];
+    a=(a^a<<12)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+    x=m[++i];
+    a=(a^a>>33)+m[i+ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+  }
+  for(i=ISAAC64_SZ/2;i<ISAAC64_SZ;i++){
+    x=m[i];
+    a=~(a^a<<21)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+    x=m[++i];
+    a=(a^a>>5)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+    x=m[++i];
+    a=(a^a<<12)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+    x=m[++i];
+    a=(a^a>>33)+m[i-ISAAC64_SZ/2]&ISAAC64_MASK;
+    m[i]=y=m[(x&ISAAC64_SZ-1<<3)>>3]+a+b&ISAAC64_MASK;
+    r[i]=b=m[y>>ISAAC64_SZ_LOG+3&ISAAC64_SZ-1]+x&ISAAC64_MASK;
+  }
+  _ctx->b=b;
+  _ctx->a=a;
+  _ctx->n=ISAAC64_SZ;
+}
+
+static void isaac64_mix(uint64_t _x[8]){
+  static const unsigned char SHIFT[8]={9,9,23,15,14,20,17,14};
+  int i;
+  for(i=0;i<8;i++){
+    _x[i]-=_x[i+4&7];
+    _x[i+5&7]^=_x[i+7&7]>>SHIFT[i];
+    _x[i+7&7]+=_x[i];
+    i++;
+    _x[i]-=_x[i+4&7];
+    _x[i+5&7]^=_x[i+7&7]<<SHIFT[i];
+    _x[i+7&7]+=_x[i];
+  }
+}
+
+
+void isaac64_init(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed){
+  _ctx->a=_ctx->b=_ctx->c=0;
+  memset(_ctx->r,0,sizeof(_ctx->r));
+  isaac64_reseed(_ctx,_seed,_nseed);
+}
+
+void isaac64_reseed(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed){
+  uint64_t *m;
+  uint64_t *r;
+  uint64_t  x[8];
+  int       i;
+  int       j;
+  m=_ctx->m;
+  r=_ctx->r;
+  if(_nseed>ISAAC64_SEED_SZ_MAX)_nseed=ISAAC64_SEED_SZ_MAX;
+  for(i=0;i<_nseed>>3;i++){
+    r[i]^=(uint64_t)_seed[i<<3|7]<<56|(uint64_t)_seed[i<<3|6]<<48|
+     (uint64_t)_seed[i<<3|5]<<40|(uint64_t)_seed[i<<3|4]<<32|
+     (uint64_t)_seed[i<<3|3]<<24|(uint64_t)_seed[i<<3|2]<<16|
+     (uint64_t)_seed[i<<3|1]<<8|_seed[i<<3];
+  }
+  _nseed-=i<<3;
+  if(_nseed>0){
+    uint64_t ri;
+    ri=_seed[i<<3];
+    for(j=1;j<_nseed;j++)ri|=(uint64_t)_seed[i<<3|j]<<(j<<3);
+    r[i++]^=ri;
+  }
+  x[0]=x[1]=x[2]=x[3]=x[4]=x[5]=x[6]=x[7]=(uint64_t)0x9E3779B97F4A7C13ULL;
+  for(i=0;i<4;i++)isaac64_mix(x);
+  for(i=0;i<ISAAC64_SZ;i+=8){
+    for(j=0;j<8;j++)x[j]+=r[i+j];
+    isaac64_mix(x);
+    memcpy(m+i,x,sizeof(x));
+  }
+  for(i=0;i<ISAAC64_SZ;i+=8){
+    for(j=0;j<8;j++)x[j]+=m[i+j];
+    isaac64_mix(x);
+    memcpy(m+i,x,sizeof(x));
+  }
+  isaac64_update(_ctx);
+}
+
+uint64_t isaac64_next_uint64(isaac64_ctx *_ctx){
+  if(!_ctx->n)isaac64_update(_ctx);
+  return _ctx->r[--_ctx->n];
+}
+
+uint64_t isaac64_next_uint(isaac64_ctx *_ctx,uint64_t _n){
+  uint64_t r;
+  uint64_t v;
+  uint64_t d;
+  do{
+    r=isaac64_next_uint64(_ctx);
+    v=r%_n;
+    d=r-v;
+  }
+  while((d+_n-1&ISAAC64_MASK)<d);
+  return v;
+}
+
+/*Returns a uniform random float.
+  The expected value is within FLT_MIN (e.g., 1E-37) of 0.5.
+  _bits: An initial set of random bits.
+  _base: This should be -(the number of bits in _bits), up to -64.
+  Return: A float uniformly distributed between 0 (inclusive) and 1
+           (exclusive).
+          The average value was measured over 2**32 samples to be
+           0.499991407275206357.*/
+static float isaac64_float_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
+  float ret;
+  int   nbits_needed;
+  while(!_bits){
+    if(_base+FLT_MANT_DIG<FLT_MIN_EXP)return 0;
+    _base-=64;
+    _bits=isaac64_next_uint64(_ctx);
+  }
+  nbits_needed=FLT_MANT_DIG-ILOGNZ_64(_bits);
+#if FLT_MANT_DIG>64
+  ret=ldexpf((float)_bits,_base);
+# if FLT_MANT_DIG>129
+  while(64-nbits_needed<0){
+# else
+  if(64-nbits_needed<0){
+# endif
+    _base-=64;
+    nbits_needed-=64;
+    ret+=ldexpf((float)isaac64_next_uint64(_ctx),_base);
+  }
+  _bits=isaac64_next_uint64(_ctx)>>64-nbits_needed;
+  ret+=ldexpf((float)_bits,_base-nbits_needed);
+#else
+  if(nbits_needed>0){
+    _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>64-nbits_needed;
+  }
+# if FLT_MANT_DIG<64
+  else _bits>>=-nbits_needed;
+# endif
+  ret=ldexpf((float)_bits,_base-nbits_needed);
+#endif
+  return ret;
+}
+
+float isaac64_next_float(isaac64_ctx *_ctx){
+  return isaac64_float_bits(_ctx,0,0);
+}
+
+float isaac64_next_signed_float(isaac64_ctx *_ctx){
+  uint64_t bits;
+  bits=isaac64_next_uint64(_ctx);
+  return (1|-((int)bits&1))*isaac64_float_bits(_ctx,bits>>1,-63);
+}
+
+/*Returns a uniform random double.
+  _bits: An initial set of random bits.
+  _base: This should be -(the number of bits in _bits), up to -64.
+  Return: A double uniformly distributed between 0 (inclusive) and 1
+           (exclusive).
+          The average value was measured over 2**32 samples to be
+           0.499990992392019273.*/
+static double isaac64_double_bits(isaac64_ctx *_ctx,uint64_t _bits,int _base){
+  double ret;
+  int    nbits_needed;
+  while(!_bits){
+    if(_base+DBL_MANT_DIG<DBL_MIN_EXP)return 0;
+    _base-=64;
+    _bits=isaac64_next_uint64(_ctx);
+  }
+  nbits_needed=DBL_MANT_DIG-ILOGNZ_64(_bits);
+#if DBL_MANT_DIG>64
+  ret=ldexp((double)_bits,_base);
+# if DBL_MANT_DIG>129
+  while(64-nbits_needed<0){
+# else
+  if(64-nbits_needed<0){
+# endif
+    _base-=64;
+    nbits_needed-=64;
+    ret+=ldexp((double)isaac64_next_uint64(_ctx),_base);
+  }
+  _bits=isaac64_next_uint64(_ctx)>>64-nbits_needed;
+  ret+=ldexp((double)_bits,_base-nbits_needed);
+#else
+  if(nbits_needed>0){
+    _bits=_bits<<nbits_needed|isaac64_next_uint64(_ctx)>>64-nbits_needed;
+  }
+# if DBL_MANT_DIG<64
+  else _bits>>=-nbits_needed;
+# endif
+  ret=ldexp((double)_bits,_base-nbits_needed);
+#endif
+  return ret;
+}
+
+double isaac64_next_double(isaac64_ctx *_ctx){
+  return isaac64_double_bits(_ctx,0,0);
+}
+
+double isaac64_next_signed_double(isaac64_ctx *_ctx){
+  uint64_t bits;
+  bits=isaac64_next_uint64(_ctx);
+  return (1|-((int)bits&1))*isaac64_double_bits(_ctx,bits>>1,-63);
+}

ccan/isaac/isaac64.h

+#if !defined(_isaac64_H)
+# define _isaac64_H (1)
+# include <stdint.h>
+
+
+
+typedef struct isaac64_ctx isaac64_ctx;
+
+
+
+#define ISAAC64_SZ_LOG      (8)
+#define ISAAC64_SZ          (1<<ISAAC64_SZ_LOG)
+#define ISAAC64_SEED_SZ_MAX (ISAAC64_SZ<<3)
+
+
+
+/*ISAAC is the most advanced of a series of pseudo-random number generators
+   designed by Robert J. Jenkins Jr. in 1996.
+  http://www.burtleburtle.net/bob/rand/isaac.html
+  This is the 64-bit version.
+  To quote:
+    ISAAC-64 generates a different sequence than ISAAC, but it uses the same
+     principles.
+    It uses 64-bit arithmetic.
+    It generates a 64-bit result every 19 instructions.
+    All cycles are at least 2**72 values, and the average cycle length is
+     2**16583.*/
+struct isaac64_ctx{
+  unsigned n;
+  uint64_t r[ISAAC64_SZ];
+  uint64_t m[ISAAC64_SZ];
+  uint64_t a;
+  uint64_t b;
+  uint64_t c;
+};
+
+
+/**
+ * isaac64_init - Initialize an instance of the ISAAC64 random number generator.
+ * @_ctx:   The ISAAC64 instance to initialize.
+ * @_seed:  The specified seed bytes.
+ *          This may be NULL if _nseed is less than or equal to zero.
+ * @_nseed: The number of bytes to use for the seed.
+ *          If this is greater than ISAAC64_SEED_SZ_MAX, the extra bytes are
+ *           ignored.
+ */
+void isaac64_init(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed);
+
+/**
+ * isaac64_reseed - Mix a new batch of entropy into the current state.
+ * To reset ISAAC64 to a known state, call isaac64_init() again instead.
+ * @_ctx:   The instance to reseed.
+ * @_seed:  The specified seed bytes.
+ *          This may be NULL if _nseed is zero.
+ * @_nseed: The number of bytes to use for the seed.
+ *          If this is greater than ISAAC64_SEED_SZ_MAX, the extra bytes are
+ *           ignored.
+ */
+void isaac64_reseed(isaac64_ctx *_ctx,const unsigned char *_seed,int _nseed);
+/**
+ * isaac64_next_uint64 - Return the next random 64-bit value.
+ * @_ctx: The ISAAC64 instance to generate the value with.
+ */
+uint64_t isaac64_next_uint64(isaac64_ctx *_ctx);
+/**
+ * isaac64_next_uint - Uniform random integer less than the given value.
+ * @_ctx: The ISAAC64 instance to generate the value with.
+ * @_n:   The upper bound on the range of numbers returned (not inclusive).
+ *        This must be greater than zero and less than 2**64.
+ *        To return integers in the full range 0...2**64-1, use
+ *         isaac64_next_uint64() instead.
+ * Return: An integer uniformly distributed between 0 and _n-1 (inclusive).
+ */
+uint64_t isaac64_next_uint(isaac64_ctx *_ctx,uint64_t _n);
+/**
+ * isaac64_next_float - Uniform random float in the range [0,1).
+ * @_ctx: The ISAAC64 instance to generate the value with.
+ * Returns a high-quality float uniformly distributed between 0 (inclusive)
+ *  and 1 (exclusive).
+ * All of the float's mantissa bits are random, e.g., the least significant bit
+ *  may still be non-zero even if the value is less than 0.5, and any
+ *  representable float in the range [0,1) has a chance to be returned, though
+ *  values very close to zero become increasingly unlikely.
+ * To generate cheaper float values that do not have these properties, use
+ *   ldexpf((float)isaac64_next_uint64(_ctx),-64);
+ */
+float isaac64_next_float(isaac64_ctx *_ctx);
+/**
+ * isaac64_next_signed_float - Uniform random float in the range (-1,1).
+ * @_ctx: The ISAAC64 instance to generate the value with.
+ * Returns a high-quality float uniformly distributed between -1 and 1
+ *  (exclusive).
+ * All of the float's mantissa bits are random, e.g., the least significant bit
+ *  may still be non-zero even if the magnitude is less than 0.5, and any
+ *  representable float in the range (-1,1) has a chance to be returned, though
+ *  values very close to zero become increasingly unlikely.
+ * To generate cheaper float values that do not have these properties, use
+ *   ldexpf((float)isaac64_next_uint64(_ctx),-63)-1;
+ *  though this returns values in the range [-1,1).
+ */
+float isaac64_next_signed_float(isaac64_ctx *_ctx);
+/**
+ * isaac64_next_double - Uniform random double in the range [0,1).
+ * @_ctx: The ISAAC64 instance to generate the value with.
+ * Returns a high-quality double uniformly distributed between 0 (inclusive)
+ *  and 1 (exclusive).
+ * All of the double's mantissa bits are random, e.g., the least significant
+ *  bit may still be non-zero even if the value is less than 0.5, and any
+ *  representable double in the range [0,1) has a chance to be returned, though
+ *  values very close to zero become increasingly unlikely.
+ * To generate cheaper double values that do not have these properties, use
+ *   ldexp((double)isaac64_next_uint64(_ctx),-64);
+ */
+double isaac64_next_double(isaac64_ctx *_ctx);
+/**
+ * isaac64_next_signed_double - Uniform random double in the range (-1,1).
+ * @_ctx: The ISAAC64 instance to generate the value with.
+ * Returns a high-quality double uniformly distributed between -1 and 1
+ *  (exclusive).
+ * All of the double's mantissa bits are random, e.g., the least significant
+ *  bit may still be non-zero even if the value is less than 0.5, and any
+ *  representable double in the range (-1,1) has a chance to be returned,
+ *  though values very close to zero become increasingly unlikely.
+ * To generate cheaper double values that do not have these properties, use
+ *   ldexp((double)isaac64_next_uint64(_ctx),-63)-1;
+ *  though this returns values in the range [-1,1).
+ */
+double isaac64_next_signed_double(isaac64_ctx *_ctx);
+
+#endif