X-Git-Url: http://git.ozlabs.org/?p=ccan;a=blobdiff_plain;f=ccan%2Fisaac%2Fisaac64.c;fp=ccan%2Fisaac%2Fisaac64.c;h=54a22095639a2730b4d9bae1fe2b6c9c1ade3e2f;hp=0000000000000000000000000000000000000000;hb=c80910455a974389d23f59e841e916427473c99e;hpb=c1bdf46aba778eb16577a23db03ef2166d9de75d diff --git a/ccan/isaac/isaac64.c b/ccan/isaac/isaac64.c new file mode 100644 index 00000000..54a22095 --- /dev/null +++ b/ccan/isaac/isaac64.c @@ -0,0 +1,250 @@ +/*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 +#include +#include +#include +#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>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>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]<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;in)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)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<>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_DIG64 + 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<>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); +}