Commit ce704f16 authored by Vladimír Čunát's avatar Vladimír Čunát

use randomness from gnutls instead of internal ISAAC

parent ebf42145
......@@ -99,7 +99,7 @@ lint:clang-scan-build:
dependencies: [] # do not download build artifacts
script:
- MAKEFLAGS="-k -j$(nproc)" SCAN_BUILD="scan-build -o scan-results --status-bugs -no-failure-reports" ./tests/clang_scan_build.sh make || true
- test "$(ls scan-results/*/report-*.html | wc -l)" = 7 # we have this many errors ATM :-)
- test "$(ls scan-results/*/report-*.html | wc -l)" = 6 # we have this many errors ATM :-)
artifacts:
when: on_failure
expire_in: '1 day'
......
......@@ -27,6 +27,7 @@ Improvements
- Dockerfile: rework, basing on Debian instead of Alpine
- policy.{FORWARD,TLS_FORWARD,STUB}: give advantage to IPv6
when choosing whom to ask, just as for iteration
- use pseudo-randomness from gnutls instead of internal ISAAC (#233)
Knot Resolver 3.1.0 (2018-11-02)
......
../../licenses/CC0
\ No newline at end of file
/**
* 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: CC0 (Public domain)
* Ccanlint:
* // We actually depend on the LGPL ilog routines, so not PD :(
* license_depends_compat FAIL
*/
#include "config.h"
#include <string.h>
#include <stdio.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/ilog\n");
return 0;
}
return 1;
}
/*Written by Timothy B. Terriberry (tterribe@xiph.org) 1999-2009.
CC0 (Public domain) - see LICENSE file for details
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"
#define ISAAC_MASK (0xFFFFFFFFU)
/* Extract ISAAC_SZ_LOG bits (starting at bit 2). */
static inline uint32_t lower_bits(uint32_t x)
{
return (x & ((ISAAC_SZ-1) << 2)) >> 2;
}
/* Extract next ISAAC_SZ_LOG bits (starting at bit ISAAC_SZ_LOG+2). */
static inline uint32_t upper_bits(uint32_t y)
{
return (y >> (ISAAC_SZ_LOG+2)) & (ISAAC_SZ-1);
}
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);
for(i=0;i<ISAAC_SZ/2;i++){
x=m[i];
a=(a^a<<13)+m[i+ISAAC_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>>6)+m[i+ISAAC_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<<2)+m[i+ISAAC_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>>16)+m[i+ISAAC_SZ/2];
m[i]=y=m[lower_bits(x)]+a+b;
r[i]=b=m[upper_bits(y)]+x;
}
for(i=ISAAC_SZ/2;i<ISAAC_SZ;i++){
x=m[i];
a=(a^a<<13)+m[i-ISAAC_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>>6)+m[i-ISAAC_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<<2)+m[i-ISAAC_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>>16)+m[i-ISAAC_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=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-ilog32_nz(_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-ilog32_nz(_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);
}
/* CC0 (Public domain) - see LICENSE file for details */
#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
/*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;