rand fails bigcrush

[andreasnoack]

This is bad news. Our random numbers should pass the tests, and a fix might not be so easy to come up with. It would be good if @ViralBShah could have a look at this.

To reproduce the fail, you'll need latest version of RNGTest and to check that the test works you'll need a genuine libRmath (not our patched version because it uses our random numbers). Basically, all p-values from the test suite should be bounded away from zero and one, but

julia> using RNGTest
julia> RNGTest.scomp_LinearComp(rand, 1, 400000, 0, 1)
(1.0,0.0418411679129036)
julia> randR()=ccall(("runif","libRmath"),Float64,(Float64,Float64),0.0,1.0)
randR (generic function with 1 method)
julia> RNGTest.scomp_LinearComp(randR, 1, 400000, 0, 1)
(0.7478822795178575,0.38169038211232587)

and it is consistent. All other tests pass for our implementation, so the problem might be dSFMT. However, I ran bigcrush in October 2012 and at that time the tests passed. I have tried to run the tests without the latest changes but it didn't help.

[ViralBShah]

Has BigCrush changed since? I believe the dSFMT authors are on github.

[StefanKarpinski]

Wait – so you think this is due to changes in dSFMT? Or our code? Or libRmath?

[andreasnoack]

I have just tried that specific test with latest version of both TestU01 and the five latest version of dSFMT and the test fails for all the versions, so I don't know what happened in the fall of 2012 when it worked. The version of TestU01 is from 2009 so I don't think that is the reason.

[andreasnoack]

Now, I think, I know what happened in 2012. At that time, we discussed the normal and exponential variates and therefore I ran BigCrush on those variates (transformed by the cdf) and not directly on rand. I have just tried to run the test that rand fails again on the normal variates and they don't fail. This could be explained by the fact that one normal actually consumes more than one uniform variate.

Because of this issue, I have tried to experiment with Random123, but I haven't managed to get nearly as fast generation as dSFMT even though they advertise that their method is faster than the Mersenne-Twister.

[ViralBShah]

Do you think that Random123 is possibly now faster than dSFMT, given that we only use genrand?

[andreasnoack]

Well. I couldn't get Random123 as fast as dSFMT when calling from Julia, but that is not the same as saying that it is not possible. As I remember it, I did a few tests in C code which suggested that Random123 is as fast as dSFMT, but the problem is that Random123 API (? I am not strong in programming terminology) is a little special. There is no source code, but only header files, so there is nothing to compile into a library. I wrote a wrapper (?), to create a library file which I could call from Julia, but that was slow. Probably that wrapper can be written better.

The code in the header files is actually quite simple, so I tried to write it in Julia. It was about the same speed as calling Random123. I considered the new vectorisation and SIMD stuff, but there are no loops in the code. It is basically control flow and XOR operations.

[ihnorton]

With a header-only library all the calls are probably inlined, but when you call into a wrapper library then the function call overhead may dominate. There isn't really a way around this right now, but it would be interesting to have a hook to "link" against an external LLVM IR file so that we could partially compile a C/C++ header with Clang and completely inline functions from external IR. Like #5046 but with just a ccall flag.

[JeffBezanson]

Any chance of doing anything about this?

[ViralBShah]

We recently disabled the array RNG in DSFMT. Would be worth rerunning this test?

[andreasnoack]

It shouldn't matter because the test only used the scalar version. I think this is property of the DSFMT, so I don't think we do much but accept the situation or change RNG. As mentioned earlier it should be noted that transforming the random variates to normality a back to uniformity removes this defect.

[ViralBShah]

Yeah, I doubt we are going to be able to do anything about it for now. Going forward, we can provide more RNGs, so that users can pick one that has the right features for the problem they are solving.

[StefanKarpinski]

I don't think providing lots of options that are all broken in different ways is an acceptable solution. The default should be a good choice that is both correct and fast. We thought DSFMT was that choice, but it's increasingly clear that it isn't. So we need something else that is actually correct and fast. If we need to make that ourselves, we can, but hopefully it doesn't come to that. In the meantime, I don't think it's appropriate to close this issue while our default RNG fails the industry standard test for good randomness.