

I was doing some testing that involved generating very large random
integers and I noticed a bug in Lua 5.3.4's math.random(). If you pass
in a range larger than 2^311 (L_RANDMAX), math.random() will generate
random integers that only have 31 bits of randomness. For example, this
code will run without errors:
for _=1,1000000 do
assert(math.random(0, 0x7fffffffffffffff) & 0xffffffff == 0)
end
print("The lowest 32 bits were always zero!")
The fix for this is pretty simple:
> ** $Id: lmathlib.c,v 1.119 2016/12/22 13:08:50 roberto Exp $
269,270c269
< luaL_argcheck(L, low >= 0  up <= LUA_MAXINTEGER + low, 1,
< "interval too large");

> luaL_argcheck(L, up  low + 1 <= L_RANDMAX, 1, "interval too large");
which makes math.random(0, 0x7fffffffffffffff) produce an
intervaltoolarge error. Also, the documentation would need to be
updated to reflect this constraint (it currently just says: "The value
nm cannot be negative and must fit in a Lua integer.").
Or, alternatively, the C code could be changed to use multiple calls to
l_random() to produce fully random integers from larger ranges, which
would introduce some additional complexity, but might be a good idea
anyway.


On Wed, Feb 28, 2018 at 07:32:47PM 0800, Bruce Hill wrote:
> I was doing some testing that involved generating very large random integers
> and I noticed a bug in Lua 5.3.4's math.random(). If you pass in a range
> larger than 2^311 (L_RANDMAX), math.random() will generate random integers
> that only have 31 bits of randomness. For example, this code will run
> without errors:
>
> for _=1,1000000 do
> assert(math.random(0, 0x7fffffffffffffff) & 0xffffffff == 0)
> end
> print("The lowest 32 bits were always zero!")
>
> The fix for this is pretty simple:
>
> > ** $Id: lmathlib.c,v 1.119 2016/12/22 13:08:50 roberto Exp $
> 269,270c269
> < luaL_argcheck(L, low >= 0  up <= LUA_MAXINTEGER + low, 1,
> < "interval too large");
> 
> > luaL_argcheck(L, up  low + 1 <= L_RANDMAX, 1, "interval too large");
up  low isn't representable when the difference is greater than
LUA_MAXINTEGER.


> On Mar 1, 2018, at 12:31 AM, William Ahern < [hidden email]> wrote:
>
>> On Wed, Feb 28, 2018 at 07:32:47PM 0800, Bruce Hill wrote:
>> I was doing some testing that involved generating very large random integers
>> and I noticed a bug in Lua 5.3.4's math.random(). If you pass in a range
>> larger than 2^311 (L_RANDMAX), math.random() will generate random integers
>> that only have 31 bits of randomness. For example, this code will run
>> without errors:
>>
>> for _=1,1000000 do
>> assert(math.random(0, 0x7fffffffffffffff) & 0xffffffff == 0)
>> end
>> print("The lowest 32 bits were always zero!")
>>
>> The fix for this is pretty simple:
>>
>>> ** $Id: lmathlib.c,v 1.119 2016/12/22 13:08:50 roberto Exp $
>> 269,270c269
>> < luaL_argcheck(L, low >= 0  up <= LUA_MAXINTEGER + low, 1,
>> < "interval too large");
>> 
>>> luaL_argcheck(L, up  low + 1 <= L_RANDMAX, 1, "interval too large");
>
> up  low isn't representable when the difference is greater than
> LUA_MAXINTEGER.
>
low >= 0 test not needed ?
Tried the assert on my Windows machine, I only get 16 bits randomness.
Probably less because it uses rand(), which is not that good.
Is there a compiler flag that let me get 32 bits randomness ?


> > > luaL_argcheck(L, up  low + 1 <= L_RANDMAX, 1, "interval too large");
>
> up  low isn't representable when the difference is greater than
> LUA_MAXINTEGER.
we have to cast all operands to LUA_UNSIGNED.
 Roberto


Sent from my iPad
> On Mar 1, 2018, at 12:31 AM, William Ahern < [hidden email]> wrote:
>
>> On Wed, Feb 28, 2018 at 07:32:47PM 0800, Bruce Hill wrote:
>> I was doing some testing that involved generating very large random integers
>> and I noticed a bug in Lua 5.3.4's math.random(). If you pass in a range
>> larger than 2^311 (L_RANDMAX), math.random() will generate random integers
>> that only have 31 bits of randomness. For example, this code will run
>> without errors:
>>
>> for _=1,1000000 do
>> assert(math.random(0, 0x7fffffffffffffff) & 0xffffffff == 0)
>> end
>> print("The lowest 32 bits were always zero!")
>>
>> The fix for this is pretty simple:
>>
>>> ** $Id: lmathlib.c,v 1.119 2016/12/22 13:08:50 roberto Exp $
>> 269,270c269
>> < luaL_argcheck(L, low >= 0  up <= LUA_MAXINTEGER + low, 1,
>> < "interval too large");
>> 
>>> luaL_argcheck(L, up  low + 1 <= L_RANDMAX, 1, "interval too large");
>
> up  low isn't representable when the difference is greater than
> LUA_MAXINTEGER.
>
I am not so sure the original version is wrong, since it does
not claim that all bits returned must be random, it just check
if the interval is too large (with lua integer). If the error message
were "not all bits are random", the patch make sense.
It might be more useful to add more random bits if the interval
is over L_RANDMAX, and not just quit . For my own use,
I patched math.random with Knuth srand64 instead of gcc rand().
rand() is not that good (at least my copy of gcc):
http://www.azillionmonkeys.com/qed/random.html


> On Mar 1, 2018, at 4:31 AM, Roberto Ierusalimschy < [hidden email]> wrote:
>
>>>> luaL_argcheck(L, up  low + 1 <= L_RANDMAX, 1, "interval too large");
>>
>> up  low isn't representable when the difference is greater than
>> LUA_MAXINTEGER.
>
> we have to cast all operands to LUA_UNSIGNED.
>
>  Roberto
Good catch and good suggestion. I believe this is correct version of what I originally proposed:
luaL_argcheck(L, (lua_Unsigned)up  (lua_Unsigned)low <= (lua_Unsigned)L_RANDMAX, 1, "interval too large");
> On Mar 1, 2018, at 4:37 AM, Albert Chan < [hidden email]> wrote:
>
> I am not so sure the original version is wrong, since it does
> not claim that all bits returned must be random, it just check
> if the interval is too large (with lua integer). If the error message
> were "not all bits are random", the patch make sense.
I think the original version is wrong because the math.random() documentation says: "... math.random returns a pseudorandom integer with uniform distribution in the range [m, n].", but the original code doesn't produce a uniform distribution in the range [m, n]. For example, math.random(2^321) only returns even numbers, which is not a uniform probability distribution over the integers in [0, 2^321] and is also very counterintuitive.
When I originally encountered this behavior, it took me a while to debug it and realize Lua was responsible for the weirdness (rather than my code), and further digging in Lua's source code to figure out why. I don't think anyone else should have to go through that process, so it would be nice to produce a helpful error message instead of silently doing something unexpected and undocumented. Either that, or have the C code support generating uniform random distributions for all integer ranges. I've written and tested a version that does this, but it's a bigger change than my original proposal, and it introduces some complexity, so I can understand if people are reluctant to make the bigger change (though I am pretty happy with this code!). Here's the change:
***************
*** 247,275 ****
static int math_random (lua_State *L) {
! lua_Integer low, up;
! double r = (double)l_rand() * (1.0 / ((double)L_RANDMAX + 1.0));
switch (lua_gettop(L)) { /* check number of arguments */
case 0: { /* no arguments */
lua_pushnumber(L, (lua_Number)r); /* Number between 0 and 1 */
return 1;
}
case 1: { /* only upper limit */
low = 1;
! up = luaL_checkinteger(L, 1);
break;
}
case 2: { /* lower and upper limits */
low = luaL_checkinteger(L, 1);
! up = luaL_checkinteger(L, 2);
break;
}
default: return luaL_error(L, "wrong number of arguments");
}
! /* random integer in the interval [low, up] */
! luaL_argcheck(L, low <= up, 1, "interval is empty");
! luaL_argcheck(L, low >= 0  up <= LUA_MAXINTEGER + low, 1,
! "interval too large");
! r *= (double)(up  low) + 1.0;
! lua_pushinteger(L, (lua_Integer)r + low);
return 1;
}
 247,295 
static int math_random (lua_State *L) {
! lua_Integer low, high;
switch (lua_gettop(L)) { /* check number of arguments */
case 0: { /* no arguments */
+ double r = (double)l_rand() * (1.0 / ((double)L_RANDMAX + 1.0));
lua_pushnumber(L, (lua_Number)r); /* Number between 0 and 1 */
return 1;
}
case 1: { /* only upper limit */
low = 1;
! high = luaL_checkinteger(L, 1);
break;
}
case 2: { /* lower and upper limits */
low = luaL_checkinteger(L, 1);
! high = luaL_checkinteger(L, 2);
break;
}
default: return luaL_error(L, "wrong number of arguments");
}
! /* random integer in the interval [low, high] */
! luaL_argcheck(L, low <= high, 1, "interval is empty");
! lua_Unsigned r, max_r; /* max_r tracks the maximum value that "r" could currently be */
! lua_Unsigned target_max_r = (lua_Unsigned)high  (lua_Unsigned)low;
! do {
! r = 0u, max_r = 0u;
! do { /* Concatenate chunks of random bits onto r until r could be >=target_max_r */
! /*
! ** Multiplying unsigned values by L_RANDMAX+1u is equivalent to left shifting by the
! ** number of random bits returned by l_rand(), since L_RANDMAX+1u is always a power of 2.
! ** (Unless L_RANDMAX+1u overflows to 0, but in that case, the loop will exit on the
! ** first iteration because max_r = ((0u * 0u)  1u) guarantees max_r >= target_max_r)
! */
! r = (r * (L_RANDMAX + 1u))  l_rand();
! max_r = (max_r * (L_RANDMAX + 1u))  L_RANDMAX;
! } while (max_r < target_max_r);
! /* If r is in exactly the range we need, it won't have bias and we can use it. */
! if (max_r == target_max_r) {
! lua_pushinteger(L, (lua_Integer)((lua_Unsigned)low + r));
! return 1;
! }
! /* Otherwise, retry if r % (target_max_r+1u) would be biased towards 0 */
! } while (r >= max_r  (max_r % (target_max_r + 1u)));
!
! r %= target_max_r + 1u; /* this will not overflow, since target_max_r < max_r */
! lua_pushinteger(L, (lua_Integer)((lua_Unsigned)low + r));
return 1;
}


On Mar 3, 2018, at 4:33 AM, Bruce Hill < [hidden email]> wrote:
> ! /* random integer in the interval [low, high] */
> ! luaL_argcheck(L, low <= high, 1, "interval is empty");
> ! lua_Unsigned r, max_r; /* max_r tracks the maximum value that "r" could currently be */
> ! lua_Unsigned target_max_r = (lua_Unsigned)high  (lua_Unsigned)low;
> ! do {
> ! r = 0u, max_r = 0u;
> ! do { /* Concatenate chunks of random bits onto r until r could be >=target_max_r */
> ! /*
> ! ** Multiplying unsigned values by L_RANDMAX+1u is equivalent to left shifting by the
> ! ** number of random bits returned by l_rand(), since L_RANDMAX+1u is always a power of 2.
> ! ** (Unless L_RANDMAX+1u overflows to 0, but in that case, the loop will exit on the
> ! ** first iteration because max_r = ((0u * 0u)  1u) guarantees max_r >= target_max_r)
> ! */
> ! r = (r * (L_RANDMAX + 1u))  l_rand();
> ! max_r = (max_r * (L_RANDMAX + 1u))  L_RANDMAX;
> ! } while (max_r < target_max_r);
> ! /* If r is in exactly the range we need, it won't have bias and we can use it. */
> ! if (max_r == target_max_r) {
> ! lua_pushinteger(L, (lua_Integer)((lua_Unsigned)low + r));
> ! return 1;
> ! }
> ! /* Otherwise, retry if r % (target_max_r+1u) would be biased towards 0 */
> ! } while (r >= max_r  (max_r % (target_max_r + 1u)));
> !
> ! r %= target_max_r + 1u; /* this will not overflow, since target_max_r < max_r */
> ! lua_pushinteger(L, (lua_Integer)((lua_Unsigned)low + r));
> return 1;
> }
Is that modified Paul Hsieh code ?
http://www.azillionmonkeys.com/qed/random.htmlDoes the do while loops guaranteed finite ?
If yes, how many loops do you expect ?
What does lua doc say about math.random with no argument ?
On my machine, math.random() only have 16 bits randomness.
(will be nice if all 53 bits are random)


Lua's math library is as good as the one that comes with your C
compiler. If your compiler is Posixcompliant, only 31 bits are random
because that is the Posix standard.
The astonishing thing is that we have knowin since 1976 how to turn a
bad random number generator into a good one, this method made it into
the second edition (1981) [1] of Donald Knuth's well known and highly
respected Art of Computer Programming, Volume 2 as Algorithm B on p32
— and most random number generators of most programming languages
still suck.
If your application requires highquality pseudorandom numbers, I
strongly urge you to implement that algorithm. Here it is in Lua.
do  closure
local X = {}
local Y
local rand = math.random
local function rand53()  integers in the range 0..(1<<53)1
return (rand(0,(1<<26)1)<<27)+rand(0,(1<<27)1)
end
function randB(init)
if init then  initialize
for k=0,255 do X[k]=rand53() end
Y = rand53()
end
j = Y>>45  highorder 8 bits
Y = X[j]
X[j] = rand53()
return Y
end
end  closure
20180303 14:24 GMT+02:00 Albert Chan < [hidden email]>:
> On Mar 3, 2018, at 4:33 AM, Bruce Hill < [hidden email]> wrote:
>
>> ! /* random integer in the interval [low, high] */
>> ! luaL_argcheck(L, low <= high, 1, "interval is empty");
>> ! lua_Unsigned r, max_r; /* max_r tracks the maximum value that "r" could currently be */
>> ! lua_Unsigned target_max_r = (lua_Unsigned)high  (lua_Unsigned)low;
>> ! do {
>> ! r = 0u, max_r = 0u;
>> ! do { /* Concatenate chunks of random bits onto r until r could be >=target_max_r */
>> ! /*
>> ! ** Multiplying unsigned values by L_RANDMAX+1u is equivalent to left shifting by the
>> ! ** number of random bits returned by l_rand(), since L_RANDMAX+1u is always a power of 2.
>> ! ** (Unless L_RANDMAX+1u overflows to 0, but in that case, the loop will exit on the
>> ! ** first iteration because max_r = ((0u * 0u)  1u) guarantees max_r >= target_max_r)
>> ! */
>> ! r = (r * (L_RANDMAX + 1u))  l_rand();
>> ! max_r = (max_r * (L_RANDMAX + 1u))  L_RANDMAX;
>> ! } while (max_r < target_max_r);
>> ! /* If r is in exactly the range we need, it won't have bias and we can use it. */
>> ! if (max_r == target_max_r) {
>> ! lua_pushinteger(L, (lua_Integer)((lua_Unsigned)low + r));
>> ! return 1;
>> ! }
>> ! /* Otherwise, retry if r % (target_max_r+1u) would be biased towards 0 */
>> ! } while (r >= max_r  (max_r % (target_max_r + 1u)));
>> !
>> ! r %= target_max_r + 1u; /* this will not overflow, since target_max_r < max_r */
>> ! lua_pushinteger(L, (lua_Integer)((lua_Unsigned)low + r));
>> return 1;
>> }
>
> Is that modified Paul Hsieh code ?
> http://www.azillionmonkeys.com/qed/random.html>
> Does the do while loops guaranteed finite ?
> If yes, how many loops do you expect ?
>
> What does lua doc say about math.random with no argument ?
>
> On my machine, math.random() only have 16 bits randomness.
> (will be nice if all 53 bits are random)
>


On Mar 3, 2018, at 8:48 AM, Dirk Laurie < [hidden email]> wrote:
> The astonishing thing is that we have knowin since 1976 how to turn a
> bad random number generator into a good one,
Don't you want the Mersenne Twister, from 1997?


On Mar 3, 2018, at 8:48 AM, Dirk Laurie < [hidden email]> wrote:
> Lua's math library is as good as the one that comes with your C
> compiler. If your compiler is Posixcompliant, only 31 bits are random
> because that is the Posix standard.
>
> The astonishing thing is that we have knowin since 1976 how to turn a
> bad random number generator into a good one, this method made it into
> the second edition (1981) [1] of Donald Knuth's well known and highly
> respected Art of Computer Programming, Volume 2 as Algorithm B on p32
> — and most random number generators of most programming languages
> still suck.
>
> If your application requires highquality pseudorandom numbers, I
> strongly urge you to implement that algorithm. Here it is in Lua.
>
> do  closure
> local X = {}
> local Y
> local rand = math.random
> local function rand53()  integers in the range 0..(1<<53)1
> return (rand(0,(1<<26)1)<<27)+rand(0,(1<<27)1)
> end
> function randB(init)
> if init then  initialize
> for k=0,255 do X[k]=rand53() end
> Y = rand53()
> end
> j = Y>>45  highorder 8 bits
> Y = X[j]
> X[j] = rand53()
> return Y
> end
> end  closure
I see some problems with above code (beside efficiency)
1. this is direct quote from Knuth book:
"shuffling methods have an inherit defect, however. They change only
the order of the generated numbers, not the numbers themselves"
2. rand53() consist of 2 rand() calls.
The 2 rand() are not independent. Even if rand() itself produce
very good random number, jamming the bits does not meant
true randomized 53 bits, and may change its period.
I tried the same method months ago by jamming 3 rand() (16 bits only)
to produce double with randomized 48 bits. But the sequence (almost)
repeat itself after few calls.
What I learn from the experience is to use published, welltested code
for math.random, and I pick a simple Knuth drand64:
do
local x = 1  random seed
function rand53()  code from Knuth AoCP 3rd edition, p 108
x = x * 6364136223846793005 + 1
return (x >> 11) * 0x1p53
end
end
I coded above in C to replace math.random()
It is much faster, and gives me 53  16 = 37 extra random bits !


We are considering xorshift128+ for Lua 5.4, so that we can have 64 bits
of randomness.
 Roberto


On Mar 3, 2018, at 10:37 AM, Roberto Ierusalimschy < [hidden email]> wrote:
We are considering xorshift128+ for Lua 5.4, so that we can have 64 bits of randomness.
I think this would be the best longterm solution, but if there's going to be a 5.3.5 bugfix release, I'd recommend including either of the patches I suggested. Also, if you're upgrading Lua's RNG, it would be nice to allow users to instantiate multiple random generators with different seeds (like Python's random.Random()). It's very useful for game development, and right now you need to use third party RNG implementations like love.math.newRandomGenerator().
Is that modified Paul Hsieh code ? http://www.azillionmonkeys.com/qed/random.htmlDoes the do while loops guaranteed finite ? If yes, how many loops do you expect ?... On my machine, math.random() only have 16 bits randomness.(will be nice if all 53 bits are random)
It's not modified Paul Hsieh code. I remembered the basic modulus/bias concept from a long time ago, and looked at some stackoverflow discussions a bit before implementing my version. The inner dowhile loop is guaranteed to run in exactly as many steps as it takes to generate enough bits for your random number. E.g. random(0, 2^161) will loop once if your platform produces 16 bits of randomness per call to l_rand(), and random(0, 2^351) will loop three times (16+16+16 = 48 >= 35). The outer dowhile loop is not guaranteed to be finite, but it has a strictly less than 50% chance of repeating the loop each time (and will typically be much less likely), so the expected number of outer loops in the worst case is <2 (and for most inputs, closer to 1).
What does lua doc say about math.random with no argument ? It says: "When called without arguments, returns a pseudorandom float with uniform distribution in the range [0,1)." I kept the existing implementation for simplicity, which uses one call to l_rand() and squashes that into [0,1). Ideally, more bits of randomness would be better, but the existing behavior here is not as much of a problem, because you normally program a bit more defensively around floating point numbers, knowing that they have limited precision and are just an approximation of Reals. If you wanted to add enough bits of randomness to fill the mantissa of an IEEE 754 floating point number, you can do something like:
lua_Unsigned max_mantissa, r = 0, max_r = 0; /* With IEEE 754, inf has 0's in all the mantissa bits and 1's everywhere else */ *(lua_Number*)&max_mantissa = 1.0/0.0; max_mantissa = ~max_mantissa; do { r = (r * (L_RANDMAX + 1u))  l_rand(); max_r = (max_r * (L_RANDMAX + 1u))  L_RANDMAX; } while (max_r < max_mantissa); lua_pushnumber(L, (lua_Number)r / ((lua_Number)max_r + 1.0)); return 1;
(Or, if you wanted to get really hacky and avoid floating point division, you could randomize the mantissa, copy the sign/exponent bits from 1.0 to get a number in [1.0,2.0), then subtract 1.0 and get a similar result to the above code)


Does the do while loops guaranteed finite ? If yes, how many loops do you expect ?
The inner dowhile loop is guaranteed to run in exactly as many steps as it takes to generate enough bits for your random number. E.g. random(0, 2^161) will loop once if your platform produces 16 bits of randomness per call to l_rand(), and random(0, 2^351) will loop three times (16+16+16 = 48 >= 35). The outer dowhile loop is not guaranteed to be finite, but it has a strictly less than 50% chance of repeating the loop each time (and will typically be much less likely), so the expected number of outer loops in the worst case is <2 (and for most inputs, closer to 1).
I would suggest you track count of rand() calls. Looking at the code, it seems to require many, many calls.
If rand() return 16 bits, inner while loop call max of 4 times. Throwing away all 4 rand calls if outer loop fail seems wasteful.
What is the chances of max_r == target_max_r ?
It just feel the code is doing too much at once. And, I am not so sure it is correct (at least about the 50% chance)
I suggest move the code as function rand_upto(), and have math_random() return low + rand_upto(high  low)
local rand, RANDMAX = math.random, 0xffff local RSIZE = RANDMAX + 1
function rand_upto(n)  return random r, 0 <= r <= n local hi = RSIZE local lo, r = n % hi, 0 if lo ~= 0 then while lo * 2 < hi do hi = hi / 2 end  minimize rand calls repeat r = rand(0, hi1) until r <= lo  remove bias end if lo == n then return r end return rand_upto((n  lo) / RSIZE) * RSIZE + r end
I tried rand_upto(1e18), rand calls average around 4 ... not bad
(1e18 = 0x0de0b6b3a7640000, so lowest 16 bits skipped rand call)


On Mar 3, 2018, at 10:01 PM, albertmcchan < [hidden email]> wrote:
> local rand, RANDMAX = math.random, 0xffff
> local RSIZE = RANDMAX + 1
>
> function rand_upto(n)  return random r, 0 <= r <= n
> local hi = RSIZE
> local lo, r = n % hi, 0
> if lo ~= 0 then
> while lo * 2 < hi do hi = hi / 2 end  minimize rand calls
> repeat r = rand(0, hi1) until r <= lo  remove bias
> end
> if lo == n then return r end
> return rand_upto((n  lo) / RSIZE) * RSIZE + r
> end
>
Just wanted to comment on above rand_upto(n), for building
an unbiased random number between 0 to n
It is meant as demo. It should be coded in C, called by math_random.
Using it as is, rand (i.e. math.random) will return result slightly biased.
I did patch my setup (luajit 1.1.8, 32bits) in C using above idea,
something unexpected happened.
I get unbiased random number(s) without performance penalty !
My guess is it avoided cost of integer to double to integer conversion.
> patched math.random now return unbiased 31bit random integer
> above rand_upto now work properly
luajit does not have 64bits integer, so above is an overkill.
I can simplify and speedup above for luajit 53bit "integer"
function rand_upto(n)  return random r = [0, n], n = [0, 0x1p53)
local lo = n % 1e9  split n, max of 2 rand calls
local r = rand(lo + 1)  1  0 <= r <= lo
if lo == n then return r end
return r + 1e9 * (rand((n  lo) / 1e9 + 1)  1)
end


Below is rand_upto(n) that work without patching any c code. It call math.random() instead of math.random(0, hi), avoiding bias.
I made an error earlier. On my machine, RAND_MAX = 0x7fff
local rand, floor = math.random, math.floor local RSIZE = 2 ^ 15  RAND_MAX + 1
function rand_upto(n)  return random r, 0 <= r <= n local hi = RSIZE local lo, r = n % hi + 1, 0 if lo > 1 then while lo + lo < hi do hi = hi / 4 end  reduce rand calls if lo >= hi then hi = hi * 2 end  scaled too far repeat r = rand() * hi until r < lo  remove bias end if lo > n then return floor(r) end return floor(r) + RSIZE * rand_upto((n  lo) / RSIZE) end

