Re: Bug in math.random() range checking

[albertmcchan <albertmcchan@...>]

On Mar 3, 2018, at 8:48 AM, Dirk Laurie <dirk.laurie@gmail.com> 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 !