From bdd10a08b179d9825c1862ff8ac94cb62799e19f Mon Sep 17 00:00:00 2001 From: Roberto Ierusalimschy Date: Mon, 26 Mar 2018 16:48:46 -0300 Subject: [PATCH] in 'random', uses high-order bits instead of low-order (better statistical properties) --- lmathlib.c | 117 +++++++++++++++++++++++++++++++---------------------- 1 file changed, 68 insertions(+), 49 deletions(-) diff --git a/lmathlib.c b/lmathlib.c index 55b44f8b..9432d405 100644 --- a/lmathlib.c +++ b/lmathlib.c @@ -1,5 +1,5 @@ /* -** $Id: lmathlib.c,v 1.126 2018/03/16 14:18:18 roberto Exp roberto $ +** $Id: lmathlib.c,v 1.127 2018/03/22 19:54:49 roberto Exp roberto $ ** Standard mathematical library ** See Copyright Notice in lua.h */ @@ -281,19 +281,19 @@ static I xorshift128plus (I *state) { /* must take care to not shift stuff by more than 63 slots */ -#define maskFIG (~(~1LLU << (FIGS - 1))) /* use FIGS bits */ -#define shiftFIG (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIG) */ +#define shiftI (64 - FIGS) /* leave FIGS bits */ +#define shiftF (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIG) */ /* ** Convert bits from a random integer into a float in the ** interval [0,1). */ static lua_Number I2d (I x) { - return (lua_Number)(x & maskFIG) * shiftFIG; + return (lua_Number)(x >> shiftI) * shiftF; } -/* convert an 'I' to a lua_Unsigned */ -#define I2UInt(x) ((lua_Unsigned)(x)) +/* convert an 'I' to a lua_Unsigned (using higher bits) */ +#define I2UInt(x) ((lua_Unsigned)((x) >> (64 - LUA_UNSIGNEDBITS))) /* convert a lua_Integer to an 'I' */ #define Int2I(x) ((I)(x)) @@ -373,40 +373,47 @@ static I xorshift128plus (I *state) { #if FIGS <= 32 -#define maskHF 0 /* do not need bits from higher half */ -#define maskLOW (~(~UONE << (FIGS - 1))) /* use FIG bits */ -#define shiftFIG (l_mathop(0.5) / (UONE << (FIGS - 1))) /* 2^(-FIG) */ +#define maskLOW 0 /* do not need bits from lower half */ +#define maskHI (~(~(lu_int32)0 >> (FIGS - 1) >> 1)) /* use FIGS bits */ +#define shiftHI 1 /* no shift */ +#define shiftF (1 / l_mathop(4294967296.0)) /* 2^(-32) */ #else /* 32 < FIGS <= 64 */ /* must take care to not shift stuff by more than 31 slots */ -/* use FIG - 32 bits from higher half */ -#define maskHF (~(~UONE << (FIGS - 33))) +/* use FIGS - 32 bits from lower half */ +#define maskLOW (~(~(lu_int32)0 >> (FIGS - 33) >> 1)) -/* use all bits from lower half */ -#define maskLOW (~(lu_int32)0) +/* use all bits from higher half */ +#define maskHI (~(lu_int32)0) -/* 2^(-FIG) == (1 / 2^33) / 2^(FIG-33) */ -#define shiftFIG ((lua_Number)(1.0 / 8589934592.0) / (UONE << (FIGS - 33))) +#define shiftHI l_mathop(4294967296.0) /* 2^32 */ + +/* 2^(-64) */ +#define shiftF ((lua_Number)(1 / (shiftHI * shiftHI))) #endif -#define twoto32 l_mathop(4294967296.0) /* 2^32 */ - static lua_Number I2d (I x) { - lua_Number h = (lua_Number)(x.h & maskHF); + lua_Number h = (lua_Number)(x.h & maskHI); lua_Number l = (lua_Number)(x.l & maskLOW); - return (h * twoto32 + l) * shiftFIG; + return (h * shiftHI + l) * shiftF; } static lua_Unsigned I2UInt (I x) { - return ((lua_Unsigned)x.h << 31 << 1) | x.l; +#if (LUA_MAXINTEGER >> 30) <= 1 +/* at most 32 bits; use only high bits */ + return ((lua_Unsigned)x.h); +#else +/* at least 33 bits */ + return ((lua_Unsigned)x.h << (LUA_UNSIGNEDBITS - 32)) | + (lua_Unsigned)x.l >> (64 - LUA_UNSIGNEDBITS); +#endif } -static I Int2I (lua_Integer n) { - lua_Unsigned un = n; - return packI((lu_int32)un, (lu_int32)(un >> 31 >> 1)); +static I Int2I (lua_Unsigned n) { + return packI((lu_int32)(n >> 31 >> 1), (lu_int32)n & ~(lu_int32)0); } #endif /* } */ @@ -420,35 +427,47 @@ typedef struct { } RanState; +/* +** Return the higher bit set in 'x' (first bit is 1). +*/ +static int higherbit (lua_Unsigned x) { + /* table of higher bits from 0 to 255 */ + static const unsigned char hb[256] = { + 0,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, + 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8, + 8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8,8 + }; + int l = 0; + while (x >= 256) { l += 8; x >>= 8; } + return l + hb[x]; +} + + /* ** Project the random integer 'ran' into the interval [0, n]. -** Because 'ran' has 2^B possible values, the projection can only be -** uniform when the size of the interval [0, n] is a power of 2 (exact -** division). To get a uniform projection into [0,lim], we first -** compute 'lim', the smallest (2^b - 1) not smaller than 'n'. If the -** random number is inside [0, n], we are done. Otherwise, we try with -** another 'ran' until we have a result inside the interval. +** To get a uniform projection into [0,n], we first compute 'shf', the +** largest number that we can right-shift 'ran' and still get numbers +** as larger as 'n'. We then shift 'ran'; if the result is inside [0, n], +** we are done. Otherwise, we try with another 'ran' until we have a +** result inside the interval. (We use right shifts to avoid the lowest +** bits of 'ran', which has poorer distributions.) */ static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n, RanState *state) { - lua_Unsigned lim = n; - if ((lim & (lim + 1)) > 0) { /* 'lim + 1' is not a power of 2? */ - /* compute the smallest (2^b - 1) not smaller than 'n' */ - lim |= (lim >> 1); - lim |= (lim >> 2); - lim |= (lim >> 4); - lim |= (lim >> 8); - lim |= (lim >> 16); -#if (LUA_MAXINTEGER >> 30 >> 2) > 0 - lim |= (lim >> 32); /* integer type has more than 32 bits */ -#endif + if (n == 0) return 0; /* special case for the unit set */ + else { + int shf = LUA_UNSIGNEDBITS - higherbit(n); + lua_assert(~(lua_Unsigned)0 >> shf >= n && /* not larger */ + ~(lua_Unsigned)0 >> shf >> 1 < n); /* largest */ + while ((ran >>= shf) > n) + ran = I2UInt(xorshift128plus(state->s)); + return ran; } - lua_assert((lim & (lim + 1)) == 0 /* 'lim + 1' is a power of 2 */ - && lim >= n /* not smaller than 'n' */ - && (lim == 0 || (lim >> 1) < n)); /* it is the smallest one */ - while ((ran & lim) > n) - ran = I2UInt(xorshift128plus(state->s)); - return ran & lim; } @@ -487,12 +506,12 @@ static int math_random (lua_State *L) { } -static void setseed (I *state, lua_Integer n) { +static void setseed (I *state, lua_Unsigned n) { int i; state[0] = Int2I(n); - state[1] = Int2I(~n); + state[1] = Int2I(0xff); /* avoid a zero state */ for (i = 0; i < 16; i++) - xorshift128plus(state); /* discard initial values */ + xorshift128plus(state); /* discard initial values to "spread" seed */ }