Re: Lua 5.4 random number generator

[Albert Chan <albertmcchan@...>]

On Mar 21, 2018, at 8:57 PM, Gé Weijers <ge@weijers.org> wrote:

I've been taking a look at the random number generator in Lua 5.4.0-work1.

I would like to suggest a different implementation of the algorithm

'project', for two reasons:

 - it does not fully eliminate the bias for small numbers that are not a

 power of two

 - more importantly: when used to flip a coin (e.g. random(0,1) =3D=3D 0)=

it

 generates a short sequence of 128 values.

The xorshift128+ RNG produces a sequence whose low-order bits are an LFSR

sequence. This is a known problem with most pseudo-random number generators

of this type. It would be better to not use the low-order bit in the way

'project' uses them.

The algorithm I use most of the time (in C99). It removes all biases, at

the cost of two divisions:

uint64_t project(uint64_t lo, uint64_t hi)

{

  assert(lo <=3D hi);

  // handle edge cases first:

  switch (hi - lo) {

  case 0:

      return lo; // one-element range

  case UINT64_MAX: // lo =3D=3D 0 and hi =3D=3D UINT64_MAX

      return xorshift128plus();

  }

  //  number of possible values

  uint64_t n =3D max - min + 1u;

  //  scale =3D=3D 2**64 / n, calculated as (2**64 - n) / n + 1u

  //  'n * scale' is the largest multiple of 'n' <=3D 2**64

  uint64_t scale =3D (uint64_t)(-n) / n + 1u;

  for (;;) {

      uint64_t r =3D xorshift128plus() / scale;

      if (r < n)

          return lo + r;

  }

}

In the case of 'n' being a power of two the high-order bits will be used

('scale' will be a power of two as well)

Similarly: the 'I2d' function should probably shift a 64-bit random value

right by 11 bits in stead of masking the value, this would get rid of the

dodgy low order bit.

BTW: the author of the xorshift128+ algorithm now recommends replacing it

with a slightly different algorithm, which is ever so slightly faster, and

has better statistical properties. It has some downsides as well, the two

low-order bits both produce an LFSR sequence, it's probably a wash.

It's here:

http://xoroshiro.di.unimi.it/

I had post the problem earlier: http://lua-users.org/lists/lua-l/2018-03/msg00332.html

Below version avoided 64-bits divisions, at the cost of possibly more generator calls.

My timings suggested this is simpler and faster, even with small n

// return unbiased random number ran, 0 <= ran <= n

// it tried ran = math.random(2 ^ (bits of n)) until ran <= n

// since the range is power of 2, it avoided the expensive 64-bits divisions

static lua_Unsigned project(lua_Unsigned ran, lua_Unsigned n, RanState *state)
{
    int bits = __builtin_clzll(n);
    while ((ran >>= bits) > n) ran = I2UInt( xorshift128plus(state->s) );
    return ran;
}