Split Mersenne-Twister states between rand() and rand(N>382)

[mbauman]

As reported on the dev mailing list, there's some sort of DSFMT state corruption after generating odd-length arrays with more than 382 elements:

julia> srand(0); rand(); x = rand(383);

julia> find(x .== x[end])
2-element Array{Int64,1}:
   2
 383

julia> find(x .== rand())
1-element Array{Int64,1}:
 3

julia> find(x .== rand())
1-element Array{Int64,1}:
 4

It looks (to my novice eyes) like rand() (which punts to an implementation in dsfmt.h:266) and rand(n) (implemented in dsfmt.c:182) end up with independent DSFMT states. I'm not sure how or why, but it only seems to get triggered after calling rand(N) with odd N > 382... which is a very specific branch in librandom.jl:120.

Edit: It seems to be much more general than what I stated above. I think rand() and rand(N) are living completely split lives with separate states; N need not be odd to trigger the bug. All that is required is that there's a call to rand() before doing the array generation (In the initial report, it was repeating calls of rand(N) with N odd — which satisfies the call to rand()).

julia> srand(0); rand()
0.8236475079774124

julia> x = rand(384);

julia> find(x .== rand())
1-element Array{Int64,1}:
 4

[mbauman]

cc: @andreasnoack

[StefanKarpinski]

This may be the worst bug we've ever had.

[JeffBezanson]

That sounds pretty extreme, but when you think about it, that really might
be true.
On Apr 18, 2014 3:38 PM, "Stefan Karpinski" [email protected]
wrote:

This may be the worst bug we've ever had.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/6573#issuecomment-40837907
.