I've described some problems with random number generation that are currently present in Nim's math module.
So I made the library "random". It was inspired by, and largely imitates the look&feel of, Python's "random" module, and even uses the same underlying implementation of Mersenne Twister (I ported it from C instead of wrapping).
So... I'm looking for feedback on the library (I am not an expert in this topic and am new to Nim).
I would be very happy to have it included with Nim.
(Updated 2015-Jan-12)
I'm not an expert here myself, so I might be wrong, but looking at mt19937ar.nim, it appears that the random integer is being converted to float by multiplications and divisions, which could lead to all sorts of weird biases.
Also, the Mersene Twister could use a lot more test cases: maybe choose a couple random seeds and copy-paste the C version's response to be compared to the nimrod version.
Another thing, although I probably sound like a broken record, but I think that modified xorshift is really cool. According to that, modified xorshift significantly outperforms Mersenne Twister on both time and quality, but not on period (although it still does very well there).
Testing is something I definitely should look into. Oh yeah, now that I think about it, testing with preselected seeds will make it deterministic, which tests always need to be.
About the underlying implementation - my argument is just going to be "Python uses it".
Although, adding other implementations can't hurt. It's just that other libraries that I found were hard to understand, with an unfriendly interface.
Note that there is already a Mersenne Twister in lib/pure/mersenne.nim. It would probably be a good idea to reuse it rather than duplicating the functionality.
Using /dev/urandom is probably not a good idea. Not only may some operating systems not support it, it may not even be available on those that do (such as when a program is running inside a chroot environment()).
Generating uniformly distributed random floating point numbers from an integer source is trickier than it sounds. See, e.g., this link for some possible issues.
genrand_int32 should probably be renamed to genrand_uint32 – the current name is confusing. I'm also not a big fan of genrand_real1 .. genrand_real3. Not only is the type called float rather than real in Nimrod, but having different functionality distinguished by an opaque integer suffix is not good API design.
Using methods is sort of pointless for your purpose. You're unlikely to operate on a heterogeneous mix of random number generators. Overloading should do what you need without incurring method dispatch overhead.
Using a closure iterator to then just call it once is even more unnecessary overhead. It also probably doesn't do what you think it does (the remaining iterations are silently discarded).
6. Yeah, that's looks like a huge problem, coming from my misunderstanding of iterators. Will look into confirming and fixing this. Aaaand fixed. Thank you.
BlaXpirit: That thing ``lib/pure/mersenne.nim`` has much, much less functionality, so you can't say I'm duplicating it.
What I'm saying is that if the existing MT implementation is lacking, then it should be improved rather than having two different MT implementations in the standard library.
BlaXpirit: In case of the pseudorandom number generator, /dev/urandom is used just for seeding (with a fallback, too).
It still introduces an unnecessary dependency on an OS feature, especially given that the fallback is lacking (epochTime() may not have a sufficient granularity, which is especially problematic when using multiple threads).
Nobody is going to use a MT for cryptography (well, if they do, they only have themselves to blame). You have to deal with basically two categories of RNG users:
Note that high-performance computers, the types of machines most likely to use RNGs for large simulations also often have a pared down OS (when you have tens of thousands of cores, you may not have the room for all the bells and whistles of a typical Linux system). Note also that a 32-bit seed is generally inadequate for scientific purposes, since repetitions are too likely with as little as ~2^16 runs (birthday paradox).
Blaxpirit: I'm not sure... Python uses it
Their implementation may or may not be correct. Contrary to popular belief, even Python programmers aren't infallible. Even so, floating point code is not something that you can always copy and paste.
Blaxpirit: Uhh no, methods do have a purpose. I don't see any other good way.
You don't need them. Procs + generics work fine. This is not a case where you need subtyping polymorphism instead of parametric polymorphism. Replace every method with a proc, remove the "abstract" methods, and use foo*[Rng](self: var Rng ...) instead of foo(self: var TRandomGenerator).
Further issues: Everything past var mersenne_twister_inst* = `` will blow up with threads. (Also, why is ``mersenne_twister_inst public – is this intentional?)
Ah, finally someone explains how to replace the usage of methods instead of just saying that methods are bad. Yeah, I'll consider using generics if they let me achieve the same result.
If mersenne_twister_inst is not public, errors appear. It doesn't like something about those generics at the end. It was intended to be not public initially. But maybe it's better this way...
And about threads... yeah... that can be bad. Not sure what to do about this.
BlaXpirit: The seed is as big as needed, if you use the function based on an array.
The point I was trying to make that even having a seed function limited to 32 bits is unlikely to be a good idea. It's useful for situations where you essentially initialize subsequent RNG states with a counter, but should in practice accommodate 64 bits.
Incidentally, the function that initializes a seed from a byte array should probably avoid the cast. You're doing a lot of extra work (and have code that is error-prone) for something that could be done more easily and cheaper with:
copyMem(addr(words[0]), addr(bytes[0]), len(bytes))
import unsigned
proc convertBytesToWords(a: openarray[uint8]): seq[uint32] =
const wsize = sizeof(uint32)
newSeq(result, (high(a) + wsize) div wsize)
for i in countup(0, high(a), 4):
let
t1 = a[i].uint32
t2 = a[i+1].uint32 shl 8
t3 = a[i+2].uint32 shl 16
t4 = a[i+3].uint32 shl 24
result[i div 4] = t1 + t2 + t3 + t4
let last = len(a) div wsize * wsize
case len(a) - last
of 1: result[result.high] = a[last].uint32
of 2: result[result.high] = a[last].uint32 + a[last+1].uint32 shl 8
of 3: result[result.high] = a[last].uint32 + a[last+1].uint32 shl 8 +
a[last+2].uint32 shl 16
else: assert len(a) == last
While this may possibly be a tiny bit slower (though it saves an allocation over your code and avoids floating point math), the rest of the seeding process is more expensive, anyway. It also has the advantage that it gives you reproducible results on both little- and big-endian machines.
I'd also in general avoid casting a seq type to a different seq type. This has two reasons:
One reason is that if Nim in the future were to get a precise reference-counting garbage collector (say, for increased interoperability with Objective C on iOS), then you'd have problems. That's not critical for personal projects, but something that probably isn't a great idea in the standard library. Casting the address of the first element to an unchecked ptr array is probably the better choice here.
The second reason is that you're breaking range checks with such a cast. This is not a problem in this particular case because you're casting a seq[uint8] to a seq[uint32], which means that range checks will not raise an error. But there are other situations where that may not be true. Again, if you need to do something like this, casting to an unchecked ptr array is probably the right thing to do.
The error that you're encountering when mersenne_twister_inst is private is probably a compiler bug. You can use:
shuffle(mersenne_twister_inst, arr)
instead of
mersenne_twister_inst.shuffle(arr)
to work around it. Usually, the rules for global symbol lookup in generic procedures are that overloaded symbols and symbols incorporated via mixin are looked up when the procedure is called in the context of the call site; other symbols are looked up when the procedure is defined in the context of the definition site. It appears that for some reason that doesn't work properly when the symbol is the left-hand side of a dot expression.
Also, while I'm looking at the code, random.random_int probably needs a revision. It is really unnecessarily expensive (floating point math, getting one byte at a time, etc.)
Here's a stab at a faster implementation (untested):
# Written in 2014 by Reimer Behrends
#
# To the extent possible under law, the author has dedicated all copyright
# and related and neighboring rights to this software to the public domain
# worldwide as per the Creative Commons CC0 Public Domain Dedication.
# This software is distributed without any warranty.
proc random_uint32(rng: var TMersenneTwister): uint32 =
rng.state.genrand_int32
proc random_uint32(rng: var TSystemRandom): uint32 =
for i in 1..sizeof(result.type):
result *= 8
result += rng.bytes_it(rng)
proc random_uint64[R](rng: var R): uint64 =
rng.random_uint32 or (rng.random_uint32 shl 32)
proc random_int*[R](rng: var R; max: Positive): Natural =
## Returns a uniformly distributed random integer ``0 <= n < max``
var mask = max.uint
for i in 0..4: # C Optimizer will unroll this
mask = mask or (mask shr (1 shl i).uint)
when sizeof(int) == 4:
while true:
result = rng.random_uint32.uint and mask
if result < max: break
else:
if max > Positive(uint32.high):
mask = mask or (mask shr 32)
while true:
result = rng.random_uint64.uint and mask;
if result < max: break
else:
while true:
result = rng.random_uint32.uint and mask
if result < max: break
A copy of the Creative Commons CC0 Public Domain Dedication can be found here.
Yes, that cast was weird, because I didn't know my way around Nimrod. But indeed, manually adding the bytes with bit shifts is much better; thanks for the heads up about endianness. Reworked.
I disagree with your remark about random_int. The perfomance can't be that bad. The one line with floating point math is necessary only because there is no operation "what position is the first non-zero bit in?" Sure, your separate implementations for 32bit and 64bit are probably faster, but what about 128bit, etc? And not all RNGs expose 32 bit numbers. I want everything to be generalized. I'm focusing on clean code more.
About mersenne_twister_inst - yes, I've seen this bug more than once. But I think I'll just leave it public anyway.
Oh wow, just switched from methods to generics. The change went extremely smoothly; didn't even get an error!
About mersenne_twister_inst - yes, I've seen this bug more than once.
It's fixed on the bigbreak branch already, I will merge it back to devel later.
BlaXPirit: I disagree with your remark about random_int. The perfomance can't be that bad.
The performance is that bad. There are several problems. One is that log2 in particular is a fairly expensive floating point operation. It makes the calculation about an order of magnitude slower. A second is that using floating point operations unnecessarily and extensively hurts hyperthreading opportunities.
With hyperthreading, each physical core supports two logical cores, who share some, but not all of the resources of the physical core, allowing them to operate in parallel if their resource usage doesn't overlap. One of the things that logical cores do share are their FP/SIMD units, so FP-intensive code automatically creates a hyperthreading bottleneck.
A third is using a loop for each byte. This is unnecessary overhead when you can do it in word-sized chunks.
BlaXpirit: Sure, your separate implementations for 32bit and 64bit are probably faster, but what about 128bit, etc?
Given that I don't know of any architecture where int (and thus Positive or Natural) would have 128 bit or more, I don't see the problem. You can, of course, generalize the code for bigger int types, but specialization for 32/64-bit still pays off.
I've pointed out the http://pcg-random.org in irc and understand your reasoning with not implementing it, but I'd like to post it here too.
<st>A bit of code review: https://github.com/BlaXpirit/nim-random/blob/master/src/random.nim#L60-L63 is buggy, certain floats will never be achieved. See http://mumble.net/~campbell/blag.txt.</st>
I'm not sure what's happening at https://github.com/BlaXpirit/nim-random/blob/master/src/random.nim#L48-L58 or how/if it compensates for non-power-of-two limits, but I have a correct (I think) implementation at https://github.com/flaviut/furry-happiness/blob/master/rand.nim#L7-L30.
Edit: oops, I didn't see the history
I've updated the library (no new release yet, so use random@#HEAD to check out the new stuff).
I've incorporated some ideas from Jehan and addressed the performance concerns.
Now it is even easier to have different RNGs under the same API. I've added Xorshift128+ and Xorshift1024*.
The main problem now is, I am extremely sloppy with testing and have no idea how to approach it. Looking for some help.
I installed random via nimble on nim 0.10.3 on linux 64 and compiled the example in the readme like so
nim c -r example.nim this works fine
nim -d:release c example.nim this fails
/.nimble/pkgs/random-0.3.0/random.nim(141, 9) Error: internal error: expr: var not init mtRandomBytes_140513 No stack traceback available