diff options
author | Roberto Ierusalimschy <roberto@inf.puc-rio.br> | 2018-03-26 16:48:46 -0300 |
---|---|---|
committer | Roberto Ierusalimschy <roberto@inf.puc-rio.br> | 2018-03-26 16:48:46 -0300 |
commit | bdd10a08b179d9825c1862ff8ac94cb62799e19f (patch) | |
tree | 562c4212e5116cdde88a6a81bbed281c40d0774d /lmathlib.c | |
parent | c5e3b2f81402773466f2afc8e87839509f65facc (diff) | |
download | lua-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.c | 117 |
1 files changed, 68 insertions, 49 deletions
@@ -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 | */ |
291 | static lua_Number I2d (I x) { | 291 | static 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 | ||
397 | static lua_Number I2d (I x) { | 398 | static 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 | ||
403 | static lua_Unsigned I2UInt (I x) { | 404 | static 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 | ||
407 | static I Int2I (lua_Integer n) { | 415 | static 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 | */ | ||
433 | static 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 | */ |
432 | static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n, | 460 | static 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 | ||
490 | static void setseed (I *state, lua_Integer n) { | 509 | static 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 | ||