aboutsummaryrefslogtreecommitdiff
path: root/lmathlib.c
diff options
context:
space:
mode:
authorRoberto Ierusalimschy <roberto@inf.puc-rio.br>2018-03-26 16:48:46 -0300
committerRoberto Ierusalimschy <roberto@inf.puc-rio.br>2018-03-26 16:48:46 -0300
commitbdd10a08b179d9825c1862ff8ac94cb62799e19f (patch)
tree562c4212e5116cdde88a6a81bbed281c40d0774d /lmathlib.c
parentc5e3b2f81402773466f2afc8e87839509f65facc (diff)
downloadlua-bdd10a08b179d9825c1862ff8ac94cb62799e19f.tar.gz
lua-bdd10a08b179d9825c1862ff8ac94cb62799e19f.tar.bz2
lua-bdd10a08b179d9825c1862ff8ac94cb62799e19f.zip
in 'random', uses high-order bits instead of low-order
(better statistical properties)
Diffstat (limited to 'lmathlib.c')
-rw-r--r--lmathlib.c117
1 files 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 @@
1/* 1/*
2** $Id: lmathlib.c,v 1.126 2018/03/16 14:18:18 roberto Exp roberto $ 2** $Id: lmathlib.c,v 1.127 2018/03/22 19:54:49 roberto Exp roberto $
3** Standard mathematical library 3** Standard mathematical library
4** See Copyright Notice in lua.h 4** See Copyright Notice in lua.h
5*/ 5*/
@@ -281,19 +281,19 @@ static I xorshift128plus (I *state) {
281 281
282/* must take care to not shift stuff by more than 63 slots */ 282/* must take care to not shift stuff by more than 63 slots */
283 283
284#define maskFIG (~(~1LLU << (FIGS - 1))) /* use FIGS bits */ 284#define shiftI (64 - FIGS) /* leave FIGS bits */
285#define shiftFIG (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIG) */ 285#define shiftF (l_mathop(0.5) / (1LLU << (FIGS - 1))) /* 2^(-FIG) */
286 286
287/* 287/*
288** Convert bits from a random integer into a float in the 288** Convert bits from a random integer into a float in the
289** interval [0,1). 289** interval [0,1).
290*/ 290*/
291static lua_Number I2d (I x) { 291static lua_Number I2d (I x) {
292 return (lua_Number)(x & maskFIG) * shiftFIG; 292 return (lua_Number)(x >> shiftI) * shiftF;
293} 293}
294 294
295/* convert an 'I' to a lua_Unsigned */ 295/* convert an 'I' to a lua_Unsigned (using higher bits) */
296#define I2UInt(x) ((lua_Unsigned)(x)) 296#define I2UInt(x) ((lua_Unsigned)((x) >> (64 - LUA_UNSIGNEDBITS)))
297 297
298/* convert a lua_Integer to an 'I' */ 298/* convert a lua_Integer to an 'I' */
299#define Int2I(x) ((I)(x)) 299#define Int2I(x) ((I)(x))
@@ -373,40 +373,47 @@ static I xorshift128plus (I *state) {
373 373
374#if FIGS <= 32 374#if FIGS <= 32
375 375
376#define maskHF 0 /* do not need bits from higher half */ 376#define maskLOW 0 /* do not need bits from lower half */
377#define maskLOW (~(~UONE << (FIGS - 1))) /* use FIG bits */ 377#define maskHI (~(~(lu_int32)0 >> (FIGS - 1) >> 1)) /* use FIGS bits */
378#define shiftFIG (l_mathop(0.5) / (UONE << (FIGS - 1))) /* 2^(-FIG) */ 378#define shiftHI 1 /* no shift */
379#define shiftF (1 / l_mathop(4294967296.0)) /* 2^(-32) */
379 380
380#else /* 32 < FIGS <= 64 */ 381#else /* 32 < FIGS <= 64 */
381 382
382/* must take care to not shift stuff by more than 31 slots */ 383/* must take care to not shift stuff by more than 31 slots */
383 384
384/* use FIG - 32 bits from higher half */ 385/* use FIGS - 32 bits from lower half */
385#define maskHF (~(~UONE << (FIGS - 33))) 386#define maskLOW (~(~(lu_int32)0 >> (FIGS - 33) >> 1))
386 387
387/* use all bits from lower half */ 388/* use all bits from higher half */
388#define maskLOW (~(lu_int32)0) 389#define maskHI (~(lu_int32)0)
389 390
390/* 2^(-FIG) == (1 / 2^33) / 2^(FIG-33) */ 391#define shiftHI l_mathop(4294967296.0) /* 2^32 */
391#define shiftFIG ((lua_Number)(1.0 / 8589934592.0) / (UONE << (FIGS - 33)))
392 392
393#endif 393/* 2^(-64) */
394#define shiftF ((lua_Number)(1 / (shiftHI * shiftHI)))
394 395
395#define twoto32 l_mathop(4294967296.0) /* 2^32 */ 396#endif
396 397
397static lua_Number I2d (I x) { 398static lua_Number I2d (I x) {
398 lua_Number h = (lua_Number)(x.h & maskHF); 399 lua_Number h = (lua_Number)(x.h & maskHI);
399 lua_Number l = (lua_Number)(x.l & maskLOW); 400 lua_Number l = (lua_Number)(x.l & maskLOW);
400 return (h * twoto32 + l) * shiftFIG; 401 return (h * shiftHI + l) * shiftF;
401} 402}
402 403
403static lua_Unsigned I2UInt (I x) { 404static lua_Unsigned I2UInt (I x) {
404 return ((lua_Unsigned)x.h << 31 << 1) | x.l; 405#if (LUA_MAXINTEGER >> 30) <= 1
406/* at most 32 bits; use only high bits */
407 return ((lua_Unsigned)x.h);
408#else
409/* at least 33 bits */
410 return ((lua_Unsigned)x.h << (LUA_UNSIGNEDBITS - 32)) |
411 (lua_Unsigned)x.l >> (64 - LUA_UNSIGNEDBITS);
412#endif
405} 413}
406 414
407static I Int2I (lua_Integer n) { 415static I Int2I (lua_Unsigned n) {
408 lua_Unsigned un = n; 416 return packI((lu_int32)(n >> 31 >> 1), (lu_int32)n & ~(lu_int32)0);
409 return packI((lu_int32)un, (lu_int32)(un >> 31 >> 1));
410} 417}
411 418
412#endif /* } */ 419#endif /* } */
@@ -421,34 +428,46 @@ typedef struct {
421 428
422 429
423/* 430/*
431** Return the higher bit set in 'x' (first bit is 1).
432*/
433static int higherbit (lua_Unsigned x) {
434 /* table of higher bits from 0 to 255 */
435 static const unsigned char hb[256] = {
436 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,
437 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,
438 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,
439 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,
440 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,
441 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,
442 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,
443 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
444 };
445 int l = 0;
446 while (x >= 256) { l += 8; x >>= 8; }
447 return l + hb[x];
448}
449
450
451/*
424** Project the random integer 'ran' into the interval [0, n]. 452** Project the random integer 'ran' into the interval [0, n].
425** Because 'ran' has 2^B possible values, the projection can only be 453** To get a uniform projection into [0,n], we first compute 'shf', the
426** uniform when the size of the interval [0, n] is a power of 2 (exact 454** largest number that we can right-shift 'ran' and still get numbers
427** division). To get a uniform projection into [0,lim], we first 455** as larger as 'n'. We then shift 'ran'; if the result is inside [0, n],
428** compute 'lim', the smallest (2^b - 1) not smaller than 'n'. If the 456** we are done. Otherwise, we try with another 'ran' until we have a
429** random number is inside [0, n], we are done. Otherwise, we try with 457** result inside the interval. (We use right shifts to avoid the lowest
430** another 'ran' until we have a result inside the interval. 458** bits of 'ran', which has poorer distributions.)
431*/ 459*/
432static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n, 460static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n,
433 RanState *state) { 461 RanState *state) {
434 lua_Unsigned lim = n; 462 if (n == 0) return 0; /* special case for the unit set */
435 if ((lim & (lim + 1)) > 0) { /* 'lim + 1' is not a power of 2? */ 463 else {
436 /* compute the smallest (2^b - 1) not smaller than 'n' */ 464 int shf = LUA_UNSIGNEDBITS - higherbit(n);
437 lim |= (lim >> 1); 465 lua_assert(~(lua_Unsigned)0 >> shf >= n && /* not larger */
438 lim |= (lim >> 2); 466 ~(lua_Unsigned)0 >> shf >> 1 < n); /* largest */
439 lim |= (lim >> 4); 467 while ((ran >>= shf) > n)
440 lim |= (lim >> 8); 468 ran = I2UInt(xorshift128plus(state->s));
441 lim |= (lim >> 16); 469 return ran;
442#if (LUA_MAXINTEGER >> 30 >> 2) > 0
443 lim |= (lim >> 32); /* integer type has more than 32 bits */
444#endif
445 } 470 }
446 lua_assert((lim & (lim + 1)) == 0 /* 'lim + 1' is a power of 2 */
447 && lim >= n /* not smaller than 'n' */
448 && (lim == 0 || (lim >> 1) < n)); /* it is the smallest one */
449 while ((ran & lim) > n)
450 ran = I2UInt(xorshift128plus(state->s));
451 return ran & lim;
452} 471}
453 472
454 473
@@ -487,12 +506,12 @@ static int math_random (lua_State *L) {
487} 506}
488 507
489 508
490static void setseed (I *state, lua_Integer n) { 509static void setseed (I *state, lua_Unsigned n) {
491 int i; 510 int i;
492 state[0] = Int2I(n); 511 state[0] = Int2I(n);
493 state[1] = Int2I(~n); 512 state[1] = Int2I(0xff); /* avoid a zero state */
494 for (i = 0; i < 16; i++) 513 for (i = 0; i < 16; i++)
495 xorshift128plus(state); /* discard initial values */ 514 xorshift128plus(state); /* discard initial values to "spread" seed */
496} 515}
497 516
498 517