3

I am trying to write a matlab function that generates a binary file containing among other things a series of 64 bit random integers. These should be of good quality, which is why I'd like to use a 64 bit mersenne twister algorithm or better. The built-in randi() function is only able to produce 32 bit numbers. I have previously generated 32 bit resultates using:

rng('shuffle', 'twister');
randi(2^32-1, 'uint32')

But this is not available in 64 bit. If I recall correctly using multiple 32 bit integers to generate a 64 bit random integer is bad practice, but if there is a good solution I'm open to it.
To make matters more difficult, I'm currently using a 32 bit windows xp machine.

Andrew Janke
  • 23,508
  • 5
  • 56
  • 85
youR.Fate
  • 796
  • 4
  • 10
  • 29
  • 1
    This is more a question about statistics, and what about combining two 32-bits numbers into one 64-bits number. – tashuhka Apr 09 '13 at 13:24
  • 1
    @youR.Fate Why is it considered bad practice to combine two 32-bit numbers into a 64-bit number? – Eitan T Apr 09 '13 at 13:32
  • as @tashuhka says, can't you just concatenate the binary representation of two 32-bit integers? I think you can prove mathematically that if the two 32-bit pseudorandoms are decent, the 64-bit pseudorandom will be decent too, no? – reverse_engineer Apr 09 '13 at 13:33
  • 1
    @Rody Oldenhis combining the 2 would be a non-issue, since i'm writing them to a binary file, so I would just write them one after the other. IIRC you lose a lot of quality if you merely concatinate 2 32 bit numbers to 1 64 bit number, since the RNG it self is flawed and you increase that further. – youR.Fate Apr 09 '13 at 13:38
  • @youR.Fate: BTW, check my updated answer; I think I've figured it out :) – Rody Oldenhuis Apr 09 '13 at 14:23
  • @youR.Fate An RNG is flawed for generating 64-bit numbers by concatenating 32-bit numbers only if the RNG was flawed at generating 32-bit numbers in the first place. – SecurityMatt Apr 10 '13 at 05:57

3 Answers3

6

It seems that support in Matlab for 64-bit random integers (and 64-bit integers in general) is still rather poor.

I suspect the "best" workaround would be to do some bithacking:

% create 32-bit uints, and cast to 64-bit
f = @() uint64( randi([0 intmax('uint32')], 'uint32') );

% bitshift and bitor to convert into a proper uint64        
R = bitor( bitshift(f(),32), f() );

or use typecast as suggested by Andrew Jake, for improved readability: 

f = @() randi([0 intmax('uint32')], 'uint32');
R = typecast([f() f()], 'uint64');

When creating more than one random number, you have to change g to: 

% bitshift and bitor: 
% ------------

% create an Nx1 uint32, and cast to 64-bit
g = @(N) uint64( randi([0 intmax('uint32')], N,1, 'uint32') );

tic
R = bitor( bitshift(g(1e7),32), g(1e7) );
toc


% typecast 
% ------------

% create a 1xN uint32, but leave the casting to typecast
g = @(N) randi([0 intmax('uint32')], 1,N, 'uint32');

tic
R = typecast([g(1e7) g(1e7)], 'uint64');   
toc

with results: 

Elapsed time is 0.717668 seconds.    % bitor/bitshift
Elapsed time is 0.705700 seconds.    % typecast w/ loop

They're equally fast, so it's really whatever you prefer.

The Mersenne Twister homepage mentions that the distribution will not change when concatenating two uint32s (thank you Andrew for the heads up), so you can indeed safely do this.

Rody Oldenhuis
  • 37,726
  • 7
  • 50
  • 96
  • If you're using a PRNG (like a PC's RNG) and the distributions look the same then, it's probably okay for most applications. +1 anyway :) – Eitan T Apr 09 '13 at 14:54
  • @EitanT: Yeah, I know that products, additions etc. of random numbers usually produce different distributions. Not sure about this case though...I *think* it's the same, no? – Rody Oldenhuis Apr 09 '13 at 15:11
  • You're adding two uniformly-distributed numbers in different ranges. If the PRNG is good enough, I guess it's probably okay. – Eitan T Apr 09 '13 at 15:16
  • 2
    The Mersenne Twister home page itself seems to say you can do this. See "Want more precision" on the FAQ page: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/efaq.html – Andrew Janke Apr 10 '13 at 03:46
  • Oh, for a terser (and IMHO more readable) type conversion, you can use `typecast` on an array of shorter ints: `f = @() randi(intmax('uint32'), 'uint32'); R = typecast( [f() f()], 'uint64' )`. – Andrew Janke Apr 10 '13 at 03:52
  • @AndrewJanke: I agree, it's more readable. However, performance might be an issue; see my update. If the aim is to create a *generic` random `uint64` generator, don't use `typedef`. Thanks for the link though! – Rody Oldenhuis Apr 10 '13 at 06:27
  • I'm using the bitor and bitshift solution. Thank you for that. Also @AndrewJanke, Thanks for the link. – youR.Fate Apr 10 '13 at 08:39
  • Typecast is plenty fast; it's a trivial operation. Anything done in a loop is going to be slower than a vectorized call, and that's a lousy loop because it's doing multiple scalar operations. The issue is that `[g(1e7) g(1e7)]` is producing an N-by-2 array, which typecast doesn't like. Switch it to `R = typecast([g(1e7); g(1e7)], 'uint64')`, or change the `g` function to return row vectors, and you'll have a vectorized implementation. Should perform the same as bitshift, because the random number generation cost should swamp everything else. – Andrew Janke Apr 10 '13 at 22:46
  • @AndrewJanke: OK, changed. The problem with the loop turned out not to be `g()` (which I didn't use anyway; see the previous edit), but `intmax` -- it's not built-in, causing the JIT accelerator to trip and leave the loop slow...Argh!! Why can't that JIT compile code recursively?! or automatically inline trivial functions!? Or do **anything** smarter than this!? – Rody Oldenhuis Apr 11 '13 at 04:27
1

Judging from this page it may be possible using the Multiplicative lagged Fibonacci generator. However I don't think it will perform better than the mersenne twister according to general standards.

Furthermore you would need a 64 bit machine to run 64 bit Matlab.

EDIT: This question would also suggest that a 64 bit version of Matlab would already produce 64 bit random integers with the default settings.

Community
  • 1
  • 1
Dennis Jaheruddin
  • 21,208
  • 8
  • 66
  • 122
  • I just tested it: `randi` 64-bit Matlab does *not* produce `uint64` values by default, and you'll get the same error message when trying (`uint64` is not supported) – Rody Oldenhuis Apr 09 '13 at 13:36
  • The 64-bit version of Matlab is just running an x64 version of the Matlab interpreter and can address more memory. It doesn't change the size of any of the Matlab datatypes or numeric functionality. But that alternate generator algorithm might work. Try doing `stream = RandStream('mlfg6331_64'); R = randi(stream, intmax('uint64'), 'uint64')`. If it's going to work, it should work in either 32-bit or 64-bit Matlab. – Andrew Janke Apr 10 '13 at 03:43
1

If you really want high quality random numbers, the 64-bit version of the Mersenne Twister is published. Here's the C implementation for it. It says it's for "64-bit machines", but I think you could compile it on a 32-bit machine, too - it's just using unsigned long long, not 64-bit pointers. Slap it in a MEX file and you could easily call it from M-code.

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c

Andrew Janke
  • 23,508
  • 5
  • 56
  • 85