diff options
author | Roberto Ierusalimschy <roberto@inf.puc-rio.br> | 2018-03-09 11:56:25 -0300 |
---|---|---|
committer | Roberto Ierusalimschy <roberto@inf.puc-rio.br> | 2018-03-09 11:56:25 -0300 |
commit | 80ae1c1c16991ef37c7360aa12da9acc68e1c0ac (patch) | |
tree | e675ffa0f22ad9cf49a1416d47e2793b5642a22e /lmathlib.c | |
parent | 40683b493474199e0d85d5dfcbb8dd6af03e6987 (diff) | |
download | lua-80ae1c1c16991ef37c7360aa12da9acc68e1c0ac.tar.gz lua-80ae1c1c16991ef37c7360aa12da9acc68e1c0ac.tar.bz2 lua-80ae1c1c16991ef37c7360aa12da9acc68e1c0ac.zip |
fairer projection of random integers into an integer interval
Diffstat (limited to 'lmathlib.c')
-rw-r--r-- | lmathlib.c | 69 |
1 files changed, 54 insertions, 15 deletions
@@ -1,5 +1,5 @@ | |||
1 | /* | 1 | /* |
2 | ** $Id: lmathlib.c,v 1.119 2016/12/22 13:08:50 roberto Exp roberto $ | 2 | ** $Id: lmathlib.c,v 1.120 2018/03/05 14:07:48 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 | */ |
@@ -273,7 +273,7 @@ static I xorshift128plus (I *state) { | |||
273 | } | 273 | } |
274 | 274 | ||
275 | 275 | ||
276 | #define mask53 (~(~0ll << 53)) | 276 | #define mask53 (~(~0LLU << 53)) |
277 | 277 | ||
278 | /* | 278 | /* |
279 | ** Convert 53 bits from a random integer into a double in the | 279 | ** Convert 53 bits from a random integer into a double in the |
@@ -283,8 +283,8 @@ static double I2d (I x) { | |||
283 | return (x & mask53) * twotomin53; | 283 | return (x & mask53) * twotomin53; |
284 | } | 284 | } |
285 | 285 | ||
286 | /* convert an 'I' to a lua_Integer */ | 286 | /* convert an 'I' to a lua_Unsigned */ |
287 | #define I2Int(x) ((lua_Integer)(x)) | 287 | #define I2UInt(x) ((lua_Unsigned)(x)) |
288 | 288 | ||
289 | /* convert a lua_Integer to an 'I' */ | 289 | /* convert a lua_Integer to an 'I' */ |
290 | #define Int2I(x) ((I)(x)) | 290 | #define Int2I(x) ((I)(x)) |
@@ -358,7 +358,7 @@ static I xorshift128plus (I *state) { | |||
358 | */ | 358 | */ |
359 | 359 | ||
360 | #define mask32 0xffffffff | 360 | #define mask32 0xffffffff |
361 | #define mask21 (~(~0 << 21)) | 361 | #define mask21 (~(~0U << 21)) |
362 | 362 | ||
363 | #define twoto32 4294967296.0 /* 2^32 */ | 363 | #define twoto32 4294967296.0 /* 2^32 */ |
364 | 364 | ||
@@ -366,12 +366,12 @@ static double I2d (I x) { | |||
366 | return ((x.x1 & mask21) * twoto32 + (x.x2 & mask32)) * twotomin53; | 366 | return ((x.x1 & mask21) * twoto32 + (x.x2 & mask32)) * twotomin53; |
367 | } | 367 | } |
368 | 368 | ||
369 | static lua_Integer I2Int (I x) { | 369 | static lua_Unsigned I2UInt (I x) { |
370 | return (((lua_Integer)x.x1 << 31) << 1) | x.x2; | 370 | return ((lua_Unsigned)x.x1 << 31 << 1) | x.x2; |
371 | } | 371 | } |
372 | 372 | ||
373 | static I Int2I (lua_Integer n) { | 373 | static I Int2I (lua_Integer n) { |
374 | return pack(n, (n >> 31) >> 1); | 374 | return pack(n, n >> 31 >> 1); |
375 | } | 375 | } |
376 | 376 | ||
377 | #endif /* } */ | 377 | #endif /* } */ |
@@ -385,9 +385,50 @@ typedef struct { | |||
385 | } RanState; | 385 | } RanState; |
386 | 386 | ||
387 | 387 | ||
388 | /* | ||
389 | ** Project the random integer 'ran' into the interval [0, n]. | ||
390 | ** Because 'ran' has 2^B possible values, the projection can only | ||
391 | ** be uniform when the size of the interval [0, n] is a power of 2 | ||
392 | ** (exact division). With the fairest possible projection (e.g., | ||
393 | ** '(ran % (n + 1))'), the maximum bias is 1 in 2^B/n. | ||
394 | ** For a "small" 'n', this bias is acceptable. (Here, we accept | ||
395 | ** a maximum bias of 0.0001%.) For a larger 'n', we first | ||
396 | ** compute 'lim', the smallest (2^b - 1) not smaller than 'n', | ||
397 | ** to get a uniform projection into [0,lim]. If the result is | ||
398 | ** inside [0, n], we are done. Otherwise, we try we another | ||
399 | ** 'ran' until we have a result inside the interval. | ||
400 | */ | ||
401 | |||
402 | #define MAXBIAS 1000000 | ||
403 | |||
404 | static lua_Unsigned project (lua_Unsigned ran, lua_Unsigned n, | ||
405 | RanState *state) { | ||
406 | if (n < LUA_MAXUNSIGNED / MAXBIAS) | ||
407 | return ran % (n + 1); | ||
408 | else { | ||
409 | /* compute the smallest (2^b - 1) not smaller than 'n' */ | ||
410 | lua_Unsigned lim = n; | ||
411 | lim |= (lim >> 1); | ||
412 | lim |= (lim >> 2); | ||
413 | lim |= (lim >> 4); | ||
414 | lim |= (lim >> 8); | ||
415 | lim |= (lim >> 16); | ||
416 | #if (LUA_MAXINTEGER >> 30 >> 2) > 0 | ||
417 | lim |= (lim >> 32); /* integer type has more than 32 bits */ | ||
418 | #endif | ||
419 | lua_assert((lim & (lim + 1)) == 0 /* 'lim + 1' is a power of 2 */ | ||
420 | && lim >= n /* not smaller than 'n' */ | ||
421 | && (lim >> 1) < n); /* it is the smallest one */ | ||
422 | while ((ran & lim) > n) | ||
423 | ran = I2UInt(xorshift128plus(state->s)); | ||
424 | return ran & lim; | ||
425 | } | ||
426 | } | ||
427 | |||
428 | |||
388 | static int math_random (lua_State *L) { | 429 | static int math_random (lua_State *L) { |
389 | lua_Integer low, up; | 430 | lua_Integer low, up; |
390 | double r; | 431 | lua_Unsigned p; |
391 | RanState *state = (RanState *)lua_touserdata(L, lua_upvalueindex(1)); | 432 | RanState *state = (RanState *)lua_touserdata(L, lua_upvalueindex(1)); |
392 | I rv = xorshift128plus(state->s); /* next pseudo-random value */ | 433 | I rv = xorshift128plus(state->s); /* next pseudo-random value */ |
393 | switch (lua_gettop(L)) { /* check number of arguments */ | 434 | switch (lua_gettop(L)) { /* check number of arguments */ |
@@ -399,7 +440,7 @@ static int math_random (lua_State *L) { | |||
399 | low = 1; | 440 | low = 1; |
400 | up = luaL_checkinteger(L, 1); | 441 | up = luaL_checkinteger(L, 1); |
401 | if (up == 0) { /* single 0 as argument? */ | 442 | if (up == 0) { /* single 0 as argument? */ |
402 | lua_pushinteger(L, I2Int(rv)); /* full random integer */ | 443 | lua_pushinteger(L, I2UInt(rv)); /* full random integer */ |
403 | return 1; | 444 | return 1; |
404 | } | 445 | } |
405 | break; | 446 | break; |
@@ -413,11 +454,9 @@ static int math_random (lua_State *L) { | |||
413 | } | 454 | } |
414 | /* random integer in the interval [low, up] */ | 455 | /* random integer in the interval [low, up] */ |
415 | luaL_argcheck(L, low <= up, 1, "interval is empty"); | 456 | luaL_argcheck(L, low <= up, 1, "interval is empty"); |
416 | luaL_argcheck(L, low >= 0 || up <= LUA_MAXINTEGER + low, 1, | 457 | /* project random integer into the interval [0, up - low] */ |
417 | "interval too large"); | 458 | p = project(I2UInt(rv), (lua_Unsigned)up - (lua_Unsigned)low, state); |
418 | r = I2d(rv); /* convert random value to a double */ | 459 | lua_pushinteger(L, p + (lua_Unsigned)low); |
419 | r *= (double)(up - low) + 1.0; /* scale it */ | ||
420 | lua_pushinteger(L, (lua_Integer)r + low); | ||
421 | return 1; | 460 | return 1; |
422 | } | 461 | } |
423 | 462 | ||