Bug in math.random() range checking

classic Classic list List threaded Threaded
16 messages Options
Reply | Threaded
Open this post in threaded view
|

Bug in math.random() range checking

Bruce Hill
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^31-1 (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
interval-too-large error. Also, the documentation would need to be
updated to reflect this constraint (it currently just says: "The value
n-m 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.

Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

William Ahern
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^31-1 (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.

Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Albert Chan

> 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^31-1 (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 ?



Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Roberto Ierusalimschy
In reply to this post by William Ahern
> > >   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

Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Albert Chan
In reply to this post by William Ahern



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^31-1 (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






Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Bruce Hill
> 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 pseudo-random 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^32-1) only returns even numbers, which is not a uniform probability distribution over the integers in [0, 2^32-1] 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;
  }




signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Albert Chan
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)

Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Dirk Laurie-2
Lua's math library is as good as the one that comes with your C
compiler. If your compiler is Posix-compliant, 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 high-quality pseudo-random 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  -- high-order 8 bits
  Y = X[j]
  X[j] = rand53()
  return Y
end
    end -- closure




2018-03-03 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)
>

Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Jay Carlson
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?

Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Albert Chan
In reply to this post by Dirk Laurie-2
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 Posix-compliant, 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 high-quality pseudo-random 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  -- high-order 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, well-tested 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) * 0x1p-53
end
end

I coded above in C to replace math.random()
It is much faster, and gives me 53 - 16 = 37 extra random bits !





Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

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

-- Roberto

Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Doug Currie
On Sat, Mar 3, 2018 at 1:37 PM, Roberto Ierusalimschy <[hidden email]> wrote:
We are considering xorshift128+ for Lua 5.4, so that we can have 64 bits
of randomness.

 Good choice! I've been using it coded in Lua...

-- e_random.lua
-- for Lua 5.3

M = {} -- the module to be returned


local function splitmix64_next (random_state)
    local x = random_state[1]
    x = x + 0x9E3779B97F4A7C15
    random_state[1] = x
    local z = x
    z = (z ~ (z >> 30)) * 0xBF58476D1CE4E5B9
    z = (z ~ (z >> 27)) * 0x94D049BB133111EB
    return z ~ (z >> 31)
end

local splitmix64_random_state = {0x5a5a5a5a5a5a5a5a5a5b, next = splitmix64_next}


local function xorshift128plus_next (random_state)
    local s1 = random_state[1]
    local s0 = random_state[2]
    local res = s0 + s1
    random_state[1] = s0
    s1 = s1 ~ (s1 << 23)
    random_state[2] = s1 ~ s0 ~ (s1 >> 18) ~ (s0 >> 5)
    return res
end

-- 

local default_random_state = {splitmix64_next(splitmix64_random_state),
                              splitmix64_next(splitmix64_random_state),
                              next = xorshift128plus_next}

local function random64 (random_state)
    if type(random_state) ~= 'table' then random_state = default_random_state end
    return random_state.next(random_state)
end


-- [0.0, 1.0]
local function random_real_53 (random_state)
    local x = random64(random_state)
    return (x >> 11) * (2.0 ^ -53)
end

-- (0.0, 1.0]
local function random_real_53ish (random_state)
    local x = random64(random_state)
    return ((x >> 11) + 1) * (2.0 ^ -53)
end

-- [0.0, 1.0]
local function random_real_63 (random_state)
    local x = random64(random_state)
    return (x >>  1) * (2.0 ^ -63)
end

M.random64    = random64
M.random_real = random_real_63

return M

-- e



Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Bruce Hill
In reply to this post by Albert Chan
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 long-term 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().

On Mar 3, 2018, at 4:24 AM, Albert Chan <[hidden email]> wrote:

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 ?
...
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 do-while 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^16-1) will loop once if your platform produces 16 bits of randomness per call to l_rand(), and random(0, 2^35-1) will loop three times (16+16+16 = 48 >= 35). The outer do-while 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 pseudo-random 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)


signature.asc (849 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Albert Chan
On Mar 3, 2018, at 6:45 PM, Bruce Hill <[hidden email]> wrote:

On Mar 3, 2018, at 4:24 AM, Albert Chan <[hidden email]> wrote:

Does the do while loops guaranteed finite ? 
If yes, how many loops do you expect ?

The inner do-while 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^16-1) will loop once if your platform produces 16 bits of randomness per call to l_rand(), and random(0, 2^35-1) will loop three times (16+16+16 = 48 >= 35). The outer do-while 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, hi-1) 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)


Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Albert Chan
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, hi-1) 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, 32-bits) 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 31-bit random integer
--> above rand_upto now work properly

luajit does not have 64-bits integer, so above is an overkill.
I can simplify and speedup above for luajit 53-bit "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
   





Reply | Threaded
Open this post in threaded view
|

Re: Bug in math.random() range checking

Albert Chan
In reply to this post by Bruce Hill
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