GitHub/lua
1
0
Fork 0
mirror of https://github.com/lua/lua.git synced 2025-02-04 06:13:04 +08:00
Code Issues Projects Releases Wiki Activity

in 'random', uses high-order bits instead of low-order

Browse Source 
(better statistical properties)
This commit is contained in:
Roberto Ierusalimschy 2018-03-26 16:48:46 -03:00
parent c5e3b2f814
commit bdd10a08b1
1 changed files with 68 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

117
lmathlib.c
View File

@ -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 ** Standard mathematical library
** See Copyright Notice in lua.h ** 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 */ /* must take care to not shift stuff by more than 63 slots */
#define maskFIG (~(~1LLU << (FIGS - 1))) /* use FIGS bits */ #define shiftI (64 - FIGS) /* leave FIGS bits */
#define shiftFIG (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIG) */ #define shiftF (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIG) */
/* /*
** Convert bits from a random integer into a float in the ** Convert bits from a random integer into a float in the
** interval [0,1). ** interval [0,1).
*/ */
static lua_Number I2d (I x) { 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 */ /* convert an 'I' to a lua_Unsigned (using higher bits) */
#define I2UInt(x) ((lua_Unsigned)(x)) #define I2UInt(x) ((lua_Unsigned)((x) >> (64 - LUA_UNSIGNEDBITS)))
/* convert a lua_Integer to an 'I' */ /* convert a lua_Integer to an 'I' */
#define Int2I(x) ((I)(x)) #define Int2I(x) ((I)(x))
@ -373,40 +373,47 @@ static I xorshift128plus (I *state) {
#if FIGS <= 32 #if FIGS <= 32
#define maskHF 0 /* do not need bits from higher half */ #define maskLOW 0 /* do not need bits from lower half */
#define maskLOW (~(~UONE << (FIGS - 1))) /* use FIG bits */ #define maskHI (~(~(lu_int32)0 >> (FIGS - 1) >> 1)) /* use FIGS bits */
#define shiftFIG (l_mathop(0.5) / (UONE << (FIGS - 1))) /* 2^(-FIG) */ #define shiftHI 1 /* no shift */
#define shiftF (1 / l_mathop(4294967296.0)) /* 2^(-32) */
#else /* 32 < FIGS <= 64 */ #else /* 32 < FIGS <= 64 */
/* must take care to not shift stuff by more than 31 slots */ /* must take care to not shift stuff by more than 31 slots */
/* use FIG - 32 bits from higher half */ /* use FIGS - 32 bits from lower half */
#define maskHF (~(~UONE << (FIGS - 33))) #define maskLOW (~(~(lu_int32)0 >> (FIGS - 33) >> 1))
/* use all bits from lower half */ /* use all bits from higher half */
#define maskLOW (~(lu_int32)0) #define maskHI (~(lu_int32)0)
/* 2^(-FIG) == (1 / 2^33) / 2^(FIG-33) */ #define shiftHI l_mathop(4294967296.0) /* 2^32 */
#define shiftFIG ((lua_Number)(1.0 / 8589934592.0) / (UONE << (FIGS - 33)))
/* 2^(-64) */
#define shiftF ((lua_Number)(1 / (shiftHI * shiftHI)))
#endif #endif
#define twoto32 l_mathop(4294967296.0) /* 2^32 */
static lua_Number I2d (I x) { 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); lua_Number l = (lua_Number)(x.l & maskLOW);
return (h * twoto32 + l) * shiftFIG; return (h * shiftHI + l) * shiftF;
} }
static lua_Unsigned I2UInt (I x) { 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) { static I Int2I (lua_Unsigned n) {
lua_Unsigned un = n; return packI((lu_int32)(n >> 31 >> 1), (lu_int32)n & ~(lu_int32)0);
return packI((lu_int32)un, (lu_int32)(un >> 31 >> 1));
} }
#endif /* } */ #endif /* } */
@ -420,35 +427,47 @@ typedef struct {
} RanState; } 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]. ** Project the random integer 'ran' into the interval [0, n].
** Because 'ran' has 2^B possible values, the projection can only be ** To get a uniform projection into [0,n], we first compute 'shf', the
** uniform when the size of the interval [0, n] is a power of 2 (exact ** largest number that we can right-shift 'ran' and still get numbers
** division). To get a uniform projection into [0,lim], we first ** as larger as 'n'. We then shift 'ran'; if the result is inside [0, n],
** compute 'lim', the smallest (2^b - 1) not smaller than 'n'. If the ** we are done. Otherwise, we try with another 'ran' until we have a
** random number is inside [0, n], we are done. Otherwise, we try with ** result inside the interval. (We use right shifts to avoid the lowest
** another 'ran' until we have a result inside the interval. ** bits of 'ran', which has poorer distributions.)
*/ */
static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n, static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n,
RanState *state) { RanState *state) {
lua_Unsigned lim = n; if (n == 0) return 0; /* special case for the unit set */
if ((lim & (lim + 1)) > 0) { /* 'lim + 1' is not a power of 2? */ else {
/* compute the smallest (2^b - 1) not smaller than 'n' */ int shf = LUA_UNSIGNEDBITS - higherbit(n);
lim |= (lim >> 1); lua_assert(~(lua_Unsigned)0 >> shf >= n && /* not larger */
lim |= (lim >> 2); ~(lua_Unsigned)0 >> shf >> 1 < n); /* largest */
lim |= (lim >> 4); while ((ran >>= shf) > n)
lim |= (lim >> 8); ran = I2UInt(xorshift128plus(state->s));
lim |= (lim >> 16); return ran;
#if (LUA_MAXINTEGER >> 30 >> 2) > 0
lim |= (lim >> 32); /* integer type has more than 32 bits */
#endif
} }
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; int i;
state[0] = Int2I(n); state[0] = Int2I(n);
state[1] = Int2I(~n); state[1] = Int2I(0xff); /* avoid a zero state */
for (i = 0; i < 16; i++) for (i = 0; i < 16; i++)
xorshift128plus(state); /* discard initial values */ xorshift128plus(state); /* discard initial values to "spread" seed */
} }