Consider using the "ziggurat" method in Random::gauss().

[Alpar Juttner]

Currently, Random::gauss() uses a Box-Muller transformation based method. Ziggurat method is reported to be a better way of generating normally distributed random numbers, though one should check if it is indeed noticeably better in practice.

[Szabó Mátyás]

This patch adds a Random::gaussZiggurat() function.

On my computer it generates much (around 5 times) faster than Box Muller. If size of int is not 32 bit, it uses less specific constants, and makes somewhat different calculations. Which in my experience makes it much slower, but still around 2 times faster than Box Muller.

Random::gauss() in this version returns Random::gaussBoxMuller() (a rename of the old method).

[Alpar Juttner]

Could you also attach the code with which the constants in random.cc was generated?

[Szabó Mátyás]

Maple sheet I used to generate the constants for random.cc

[Szabó Mátyás]

I hope it's understandable. :)

[Alpar Juttner]

3ae5695818c3-b33813f266e3.patch​ contains two changesets:

• [3ae5695818c3] is a somewhat reworked version of your patch
  • Put Ziggurat related constants into a separate file (random_ziggurat.cc)
  • Broke those array to one number per lines
  • Fix the indentation of random.h
• [b33813f266e3] contains some more doc improvements and formatting.

[Alpar Juttner]

I have two questions to be answered before merging to the main:

• I don't like the if(sizeof(int) == 4) {...} condition. Isn't there any better solution for that? E.g. check this at config time. As I recall, similar distinction has to be made inside the Mersenne-Twister algorithm itself. Balazs, how do you handle this?
• Shall we use Ziggurat by default instead of Box-Muller?
  • Pros (a single but very strong one)
    • Ziggurat is much faster than Box-Muller
  • Cons
    • rnd.gauss() will generate a different sequence of numbers using the new releases. This is a very slight backward incompatibility.
    • It increases the executable size with some 15kB.
    • Ziggurat method accesses the above memory block at random places, which might worsen the cache performance when its call is embedded into another algorithm. (This is a mere speculation).

[Szabó Mátyás]

The problem that arises here, (that is: whether size or speed has priority) might occur in other algorithms as well. Maybe (as a future user) it would be useful to have an overall option to chose which is more important, like "#define PRIORITY_SPEED" or something... Just an idea :)

[Alpar Juttner]

By the way, what does if(sizeof(int) == 4) mean? Is there any contemporary system where int is not 4 byte? If yes, how big it is?

[Szabó Mátyás]

Replying to alpar:

By the way, what does if(sizeof(int) == 4) mean? Is there any contemporary system where int is not 4 byte? If yes, how big it is?

I know of none. However C++ standard enables it, and the algorithm needs the int to be 32 bits. (For the use of the previously calculated constants.)

[Alpar Juttner]

Sorry, but I'm still confused.

As far as I understand, we have a special (faster) implementation when sizeof(int) == 4. But what does sizeof(int) != 4 mean? E.g. will the algorithm work properly with 16 bit integers? Or does it (basically) refer to 64 bit architectures?

[Szabó Mátyás]

Replying to alpar:

Sorry, but I'm still confused.

As far as I understand, we have a special (faster) implementation when sizeof(int) == 4. But what does sizeof(int) != 4 mean? E.g. will the algorithm work properly with 16 bit integers? Or does it (basically) refer to 64 bit architectures?

It will. Or at least it should (I didn't test it). Should work with everything greater than or equal to 8 bit.

The only thing we use about 32 bit integers, is that the biggest absolute value possible is 2^{31. So we can multiply our consts by 2}31 (and by 2^{(-31)), and use them instead. Which makes it faster somehow...}

[Szabó Mátyás]

Sorry about the previous one. Forgot to preview first.

Nice ^feature by the way. :)

Replying to alpar:

Sorry, but I'm still confused.

As far as I understand, we have a special (faster) implementation when sizeof(int) == 4. But what does sizeof(int) != 4 mean? E.g. will the algorithm work properly with 16 bit integers? Or does it (basically) refer to 64 bit architectures?

It will. Or at least it should (I didn't test it). Should work with everything greater than or equal to 8 bit.

The only thing we use about 32 bit integers, is that the biggest absolute value possible is 2³¹. So we can multiply our consts by 2³¹ (and by 2^(-31)), and use them instead. Which makes it faster somehow...

[Smartmil8]

​about.me ​purevolume.com

[Smartmil8]

​wardsauto.com ​winsupersite.com