17

The C++ standard library does contain a set of random distribution classes, among other things a Mersenne Twister engine.

What I'd like to do is to generate the same random number distribution, using the same seed for C++ and Mathematica.

For instance, generate 10 random numbers using the mersenne twister engine (C++):

std::vector<size_t> a1( 10 );
std::mt19937 gen(42);
std::generate(a1.begin(), a1.end(), gen);

/* {1608637542, 3421126067, 4083286876, 787846414, 3143890026, 3348747335, 
    2571218620, 2563451924, 670094950, 1914837113} */

and now I generate 10 random integers in Mathematica using the MersenneTwister engine:

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"]
RandomInteger[{0, 2^32 - 1}, 10]]

(* {3012359023, 1649004928, 188383571, 1935467488, 3494372723, 668236121, 
    1292572136, 98984411, 2487091843, 3121826951} *)

Apparently not the same. The mersenne twister in the C++ standard library uses the following default settings to generate the distribution:

enter image description here

The question is now: How much influence do I have on Mathematica's random number generator, so that I can overwrite some specific settings to generate the same set of numbers for a given seed?

EDIT: MathLink implementation MT19937Range:

The .tm file:

:Begin:
:Function:      mersenneRange
:Pattern:       MT19937Range[seed_, len_]
:Arguments:     { seed, len }
:ArgumentTypes: { Integer, LongInteger }
:ReturnType:    Manual
:End:

:Evaluate:      MT19937Range::usage = ""

The .cpp file:

#include <random>

using namespace std;

void mersenneRange( int seed, long length )
{
    mt19937 gen(seed);

    MLNewPacket(stdlink);
    MLPutFunction(stdlink, "List", length);    
    for(long x = 1; x <= length; x++)
        MLPutLongInteger(stdlink, gen());
}

MLPutFunction expects as a third argument a int. So implicit conversion here with loss of precision. Is there a better alternative?

The Mathematica side:

link = Install["wherever the binary is...",  LinkMode -> Launch]

MT19937Range[42, 10]

(* {1608637542, 3421126067, 4083286876, 787846414, 3143890026, 3348747335, 2571218620,
   2563451924, 670094950, 1914837113} *)

Uninstall[link]
Stefan
  • 5,347
  • 25
  • 32
  • Yes @OleksandrR my fault. Was a typo... – Stefan Feb 13 '13 at 15:16
  • You can define your own random number generation method for Mathematica, see http://tpfto.wordpress.com/2012/02/12/on-emulating-the-texas-instruments-random-number-generator/ for an example. The docs describe the API to do this. You can write a LibraryLink function which calls your C++'s RNG. – Szabolcs Feb 13 '13 at 19:15
  • See here for the three different Mersenne twister implementations Mathematica includes, and also for defining your own custom RNG method. – Szabolcs Feb 13 '13 at 19:18
  • 2
    @Szabolcs thank you for the links. At the moment I'm writing on a skeleton implementation for integration into the RNG framework. But, nevertheless, this wasn't the question. The mt19937 is well defined and I just want to know what parameter values the folks at wolfram are using and if I can change them...like option values. The values I've listed are just parameter values for the typedef in the C++ standard library. – Stefan Feb 13 '13 at 19:38
  • 1
    @Stefan Are you sure it wouldn't also depend on hardware parameters, e.g. the size of a machine integer? You might not have control over this in Mathematica, so just using the C++ generator may be better ... I think the most complete docs on RNG methods is what I linked. It doesn't describe suboptions for MersenneTwister there, but there may be some undocumented ones. I'd just write to support and ask about it. They may be slow to reply, but they're pretty good at providing this type of information. If there are suboptions, they'll let you know. – Szabolcs Feb 13 '13 at 19:42
  • Re your edit "Is there a better alternative?" --> The better alternative is using LibraryLink. This gives you direct access to the kernel's memory and the data transfer overhead will be minimal. With MathLink, the data transfer overhead is so large that you really don't want to send more elements that what fit in a 32-bit int. Whether it is possible to put more than $2^{31}$ elements through MathLink, I do not know. – Szabolcs Feb 15 '13 at 15:34
  • @Szabolcs I see. Another alternative would be to use a buffered transfer for MathLink. So the general advise is to use LibraryLink with Mma >= 8 and forget MathLink? – Stefan Feb 15 '13 at 15:43
  • Actually it depends on the situation. Also here I must state that while I have used both LibraryLink and MathLink, I don't have extensive experience with either. The big advantage of LibraryLink is that there's no data transfer overhead. You access kernel memory directly. It's also easier to learn and to set up. The disadvantage is that your functions are loaded directly into the kernel. If they crash, the kernel crashes. It's also less flexible: you can only work with numbers, but not symbols (to be precise you can use the MathLink API in LibraryLink functions too). – Szabolcs Feb 15 '13 at 15:48
  • understood. Thank you. May I ask you if you know how to return a True/False? I'd like to check something out and need this return type for Select[] to pick the elements... – Stefan Feb 15 '13 at 15:59
  • If I get the time I'll show you how to convert your MathLink program to LibraryLink, but right now I can't do that. Maybe tonight. – Szabolcs Feb 15 '13 at 16:12
  • I did not yet get around to doing this. Did you manage to get it working with LibraryLink? – Szabolcs Feb 19 '13 at 22:43
  • Hey @Scabolcs. I did not either. So far I'm well with MathLink for the moment. LibraryLink looks to me really weird. Due to the Fibonacci range product question, where I was hacking around using C++11 techniques to calculate the results I gave JLink a try and implemented there a full blown version with arbitrary precision support plus parallel Karatsuba, ToomCook3 and SchönhageStrassen. That was quite fun, but I do not dare to post this, since I obv. tend to exceed the answer limits ;) So I wait for the question on how to do this in Java/JLink. Cheers – Stefan Feb 22 '13 at 08:57
  • And since yesterday I do have a working implementation for Takahashi Fibonacci algorithm. – Stefan Feb 22 '13 at 08:58

1 Answers1

11

If you need to guarantee that Mathematica will generate the same random numbers as a particular C++ implementation, then it's easier to just use the C++ generator in Mathematica as well.

Mathematica allows defining your own random number generation method, that can be used in any of Mathematica's random functions.

It should be possible, and relatively quick, to write a LibraryLink function that uses C++'s RNG, and then use that to define a custom RNG method in Mathematica. See here on getting started with LibraryLink. It's also possible to put the short C++ program in a string and avoid explicitly creating a separate file.

Community
  • 1
Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
  • Interesting. LibraryLink was news to me, but using the same parameter sets plus the same seed hast to provide the same answers. On the standard library side these values and the imementation are open. I'm not interested how, at wolfram, the guys did implement this. I was asking because I'd like to verify some algorithms I wrote in C++, in Mathematica to proof they work correct. They use a different parameter set. Why? Can I change it as a user? – Stefan Feb 13 '13 at 19:49
  • @Stefan they could even be using the same parameters--if the method for producing particular values from the RNG output differs, the results will still be different. I would think this is a much more likely explanation than that they decided to use the Mersenne twister in an unusual way. – Oleksandr R. Feb 13 '13 at 19:56
  • @OleksandrR. hmm...MT19937 is a pseudo-random generator. You enter a specific seed value and it is guaranteed that the distribution will always be identical...if not, well, flawed. I don't say, they've used it in a unusual way. Maybe it is not a MT19937 at all, which would make me sad. Btw. I've already exported the distribution to a csv file days before, in order to compare the results of my algorithms with those from Mathematica to ensure I'm correct. Don't you think that you should be allowed to have influence on parameters within the RNG framework? Well, I think so. – Stefan Feb 13 '13 at 21:09
  • @Stefan To put it another way: if you need a quick solution right now, I suggest you use C++'s RNG as I described. If you really want to know if you can set parameters, the best way to find out is asking support@wolfram.com. This will take some time though and there may not be such parameters at all. Or if there are, you may not have access to all of them, e.g. word length. If they reply you, please post the reply here as an answer. – Szabolcs Feb 13 '13 at 21:12
  • 3
    @Stefan you missed my point. MT19937 generates pseudorandom strings of 32 bits at a time. How that output is used by Mathematica to produce a result drawn from a particular distribution is completely unspecified and could easily account for the difference. I highly doubt that it is giving you the raw output from the underlying RNG. – Oleksandr R. Feb 14 '13 at 04:20