add file for the isaac64 random number generator
This is the version from the Comprehensive C Archive Network, licensed under the CC0 "No Rights Reserved" Creative Common License. http://ccodearchive.net/info/isaac.html
This commit is contained in:
Patric Mueller
131
include/isaac64.h
Normal file
131
include/isaac64.h
Normal file
@@ -0,0 +1,131 @@
|
||||
/* CC0 (Public domain) - see LICENSE file for details */
|
||||
#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
|
||||
255
src/isaac64.c
Normal file
255
src/isaac64.c
Normal file
@@ -0,0 +1,255 @@
|
||||
/*Written by Timothy B. Terriberry (tterribe@xiph.org) 1999-2009
|
||||
CC0 (Public domain) - see LICENSE file for details
|
||||
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"
|
||||
|
||||
|
||||
#define ISAAC64_MASK ((uint64_t)0xFFFFFFFFFFFFFFFFULL)
|
||||
|
||||
/* Extract ISAAC64_SZ_LOG bits (starting at bit 3). */
|
||||
static inline uint32_t lower_bits(uint64_t x)
|
||||
{
|
||||
return (x & ((ISAAC64_SZ-1) << 3)) >>3;
|
||||
}
|
||||
|
||||
/* Extract next ISAAC64_SZ_LOG bits (starting at bit ISAAC64_SZ_LOG+2). */
|
||||
static inline uint32_t upper_bits(uint32_t y)
|
||||
{
|
||||
return (y >> (ISAAC64_SZ_LOG+3)) & (ISAAC64_SZ-1);
|
||||
}
|
||||
|
||||
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);
|
||||
for(i=0;i<ISAAC64_SZ/2;i++){
|
||||
x=m[i];
|
||||
a=~(a^a<<21)+m[i+ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
x=m[++i];
|
||||
a=(a^a>>5)+m[i+ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
x=m[++i];
|
||||
a=(a^a<<12)+m[i+ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
x=m[++i];
|
||||
a=(a^a>>33)+m[i+ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
}
|
||||
for(i=ISAAC64_SZ/2;i<ISAAC64_SZ;i++){
|
||||
x=m[i];
|
||||
a=~(a^a<<21)+m[i-ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
x=m[++i];
|
||||
a=(a^a>>5)+m[i-ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
x=m[++i];
|
||||
a=(a^a<<12)+m[i-ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
x=m[++i];
|
||||
a=(a^a>>33)+m[i-ISAAC64_SZ/2];
|
||||
m[i]=y=m[lower_bits(x)]+a+b;
|
||||
r[i]=b=m[upper_bits(y)]+x;
|
||||
}
|
||||
_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-ilog64_nz(_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-ilog64_nz(_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);
|
||||
}
|
||||
Reference in New Issue
Block a user