What function generates random numbers in a compiled function?

12

I am writing a simulation for 2D random walk. I wrote several version of code, and found that the speed is like this (fastest to slowest): C++, Mathematica compiled procedural code, Mathematica functional code. However, the rand() in C is not good enough for the experiment, and it is said that Mathematica uses cellular automation to generate random numbers, which have better distribution, so I wish to use the latter one.

There are two cases, the first one is just plainCompile, the second one uses the option: CompilationTarget -> "C", which is about 4 times faster than the first one. But the documentation states something scary (for the second case): the function is converted to C, compiled to DLL and linked back. I am wondering that if the compiled program uses the wolfram random algorithm or the rand() provided by C in those cases.

PS: I am using RandomInteger and RandomReal for random number generation.

3Why don't you use a C++ library, e.g. Boost.Random? It is not generally true that rand() from the C standard library is not good enough. Some implementations are good, some are not so good. The problem is that the implementation you get depends on the system/compiler you use. If you use Boost.Random or a similar library you'll know which precise implementation you have and you'll even have a choice of several different implementations. – Szabolcs – 2015-07-02T11:58:51.310

2That said, your question is valid and interesting, but given your motivation for asking it: I'd just use a C++ library. – Szabolcs – 2015-07-02T11:59:47.463

3What I sometimes do is generate a huge block of random number in Mathematica and send this (as a block) to a compiled Mathematica routine. This is quite often much faster than generating the random numbers one by one in the compiled routine. – Sjoerd C. de Vries – 2015-07-02T12:01:49.270

@Szabolcs The actual reason I wanted to switch is that I need to use the data generated by C++, and plot or analysis it in Mathematica, I have to let the program output a file and let mathematica imports it, it is more complicated in Mathematica. – vapor – 2015-07-02T12:02:49.600

1@Felix You can make C++ functions accessible from Mathematica through LibraryLink. This does take a bit of work, but once you've got used to it it's not so bad. This is now the usual way for me to call C++ code (I hate writing to files using newly invented formats and importing the data back). Of course if it's fast enough in pure Mathematica that's much less trouble. – Szabolcs – 2015-07-02T12:06:11.837

2

Or even use C++11 <random>? http://www.cplusplus.com/reference/random/

– dr.blochwave – 2015-07-02T12:24:36.733

@SjoerdC.deVries That's a good idea. May I know more details on passing that block of numbers to the function? I think extracting an element from a huge array may take sometime. – vapor – 2015-07-02T12:36:43.097

I think it is unlikely that whatever Mathematica implements will be dramatically better or worse (by whatever metric) than Mersenne twister, and I certainly would not rely on some rumored superiority for the correctness of my results! If it is very important, you should check the code for the RNG yourself. Thus, it comes back to @Szabolcs's suggestion: use a suitable C++ RNG library that you know has passed the Diehard and Bigcrush tests. See also: Quality of random numbers

– Oleksandr R. – 2015-07-02T13:00:16.003

@Felix Extracting an element from an array, i.e. numbers stored sequentially in memory in a flat format, takes the same amount of time regardless of the size of the array. – Szabolcs – 2015-07-02T13:02:37.650

1You may also like to note that many modern CPUs have a built-in physical random number generator. It may not be suitable (by itself) for cryptographic purposes because its implementation is not totally open and robust against hijacking, but since it is based on thermal noise, there are good physical reasons to suppose that it will produce good-quality random output. – Oleksandr R. – 2015-07-02T13:13:41.183

14

This is not a full answer, but too long for a comment.

We can try to look at the C code that Mathematica generates.

Needs["CCodeGenerator"]

cf = Compile[{{n, _Integer}}, RandomReal[1, n]]
CCodeStringGenerate[cf, "fun"]


Excerpt from the output:

...
FP0 = funStructCompile->getFunctionCallPointer("RandomReals");
...
err = FP0(libData, 3, FPA, FPA[3]);/*  RandomReals  */
...


This is the function that's being called. What this is, we can only guess, but I'd be very surprised if it were simply using the C standard library's random number generator. It very likely uses a cross platform implementation contained within the Mathematica kernel.

Update

Here's proof that the same random number generator is being used regardless of whether the function is compiled and regardless of the compilation target.

cf =
Compile[{{n, _Integer}}, SeedRandom[123]; RandomReal[1, n]];

cf2 =
Compile[{{n, _Integer}}, SeedRandom[123]; RandomReal[1, n],
CompilationTarget -> "C"];

cf[5]
(* {0.45571914, 0.97782592, 0.94321464, 0.96221575, 0.3023479} *)

cf2[5]
(* {0.45571914, 0.97782592, 0.94321464, 0.96221575, 0.3023479} *)

SeedRandom[123]; RandomReal[1, 5]
(* {0.45571914, 0.97782592, 0.94321464, 0.96221575, 0.3023479} *)


I haven't tried other random number generation methods than the default. You might want to experiment with other options as well (the Method option of SeedRandom).

Thanks! And this proves what Sjoerd C. de Vries said: the uncompiled version is faster. – vapor – 2015-07-02T12:40:54.063

Seems to work with different method of random number generation... – dr.blochwave – 2015-07-02T12:58:42.857

2Try this with something other than a uniform distribution, e.g. cf = Compile[{{n, _Integer}}, RandomVariate[ChiSquareDistribution[5]]]` and you see it calls back to Mathematica, – Sjoerd C. de Vries – 2015-07-02T13:27:19.673

2

– dr.blochwave – 2015-07-02T14:16:47.713