Toby Ewing wrote:
Frank Heckenbach wrote:
Toby Ewing wrote:
I've been working with the Mersenne Twister random number generator for a while. I translated it from the original C to Pascal several years ago, and the GPC team incorporated it into the GPC unit. But since then the original authors have found and fixed a bug in their initilization procedures, and other users have suggested various speedups.
I've now updated my original translation, validated it against the authors' test output, performed some other minor tests, and done what small optimizations I know how to do. Source for the unit and a test program are appended below.
I've briefly looked at your code, and it seems there are some off-by-one errors, e.g.:
for kk := mtN-mtM to mtN-1 do begin
(where the C version says `<N-1'). I think there are more.
There are indeed several such differences, but they're not errors as such. I used the Pascal [1..n] convention, rather than staying with the C [0..n-1] convention.
Then these lines seem to be wrong:
type MTary = array[0..mtN] of MTcard; magary = array[0..1] of MTcard;
for kk := 0 to mtN-mtM-1 do begin
(perhaps more)
Well, `[0 .. mtN]' is one too large for C's `[mtN]' in any case. Maybe the extra element is just declared and computed but never used, so the results are the same, but it's still a bit strange.
Also, if you'd like to have this code included with GPC, it would be good to use the same coding style and interfaces (including function names such as `Default_RandInt' which are called through pointers in the RTS). I don't have the time to do the reformatting myself now.
The basics - initialization from a given random number seed, and the "rand" that returns a real in [0..1) - are currently part of the gpc.pas unit.
The actual code is in rts/random.pas in the GPC sources. gpc.pas just contains the interface.
I'd be happy to have the current code included - that was my purpose in submitting it - but I don't understand the coding and interface conventions you're using. If you can point me to some clear guidelines, I could give it a shot.
See the info file GPCS that comes with GPC (long) and/or just follow the style and interfaces in random.pas. I'd also appreciate if you keep the identifier names as far as possible so it's easy to see what's actually changed.
If you want to add more interfaces, I'd agree if they're nontrivial enough -- but in contrast e.g. `randrEx1 (range)' is equivalent to (and in fact does the same as) `range * Random' when `Random' calls `Default_RandReal' internally, so I don't think it justifies another RTS function.
I also doubt whether the bit-fiddling in `randrInt' is really better than our `Default_RandInt' code (which looks more complicated, but only because it also provides a 64 bit range -- though it may be worthwhile to add a "fast path" for MaxValue <= RandomRange there). Assuming that the random operation is the most costly part, using `mod' seems better than `and'ing out bits -- how much, of course, depends on the number, e.g. for 0 .. 32, the `and' version will require almost two calls on average (64/33) whereas the `mod' version takes just slightly more than one (1 + 4/(2^32-4)).
At a quick glance, also `randNorm' seems problematic since if `bsecond' is set, it ignores its arguments.
The code that I started with, the C++ abyss from which I rescued some civilized Pascal, claims to be faster than what I wrote. It uses the C and C++ directives "inline" and "register", and some arcane C constructions that are better avoided in clear code.
Yes, GPC provides them as attributes, but with `-O3' for `inline' and `-O' for `register', it also does it automatically, so they're recommended only in very special cases.
The other changes are
- I split the basic number generator into two parts - one that is executed for
each number, the other of which is executed only once every 624 calls. In my tests, this resulted in a 20+% speedup, presumably because the memory overhead for the more frequent function call is greatly reduced.
Yes, this may be worthwhile (with inlining).
Frank