COIN-OR::LEMON - Graph Library

Opened 9 years ago

Last modified 7 years ago

#105 new enhancement

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

Reported by: alpar Owned by: alpar
Priority: major Milestone:
Component: core Version: hg main
Keywords: random Cc: deba
Revision id:

Description

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.

Attachments (5)

2fefbb0c8a28.patch (58.2 KB) - added by SzOtYeSz 7 years ago.
Ziggurat.mw (113.8 KB) - added by SzOtYeSz 7 years ago.
Maple sheet I used to generate the constants for random.cc
3ae5695818c3-b33813f266e3.patch (71.0 KB) - added by alpar 7 years ago.
Virtual Campus Supercomputing Center.png (14.0 KB) - added by Slavon 3 years ago.
about.me purevolume.com
developedBy.jpg (4.0 KB) - added by Slavon 3 years ago.
wardsauto.com winsupersite.com

Download all attachments as: .zip

Change History (17)

comment:1 Changed 9 years ago by alpar

  • Milestone LEMON 1.1 release deleted

Changed 7 years ago by SzOtYeSz

comment:2 Changed 7 years ago by SzOtYeSz

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).

comment:3 Changed 7 years ago by alpar

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

Changed 7 years ago by SzOtYeSz

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

comment:4 Changed 7 years ago by SzOtYeSz

I hope it's understandable. :)

Changed 7 years ago by alpar

comment:5 Changed 7 years ago by alpar

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.

comment:6 follow-up: Changed 7 years ago by alpar

  • Cc deba added
  • Priority changed from minor to major

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).

comment:7 in reply to: ↑ 6 Changed 7 years ago by SzOtYeSz

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 :)

comment:8 follow-up: Changed 7 years ago by 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?

comment:9 in reply to: ↑ 8 Changed 7 years ago by SzOtYeSz

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.)

comment:10 follow-ups: Changed 7 years ago by 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?

comment:11 in reply to: ↑ 10 Changed 7 years ago by SzOtYeSz

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 231. So we can multiply our consts by 231 (and by 2(-31)), and use them instead. Which makes it faster somehow...

comment:12 in reply to: ↑ 10 Changed 7 years ago by SzOtYeSz

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 231. So we can multiply our consts by 231 (and by 2(-31)), and use them instead. Which makes it faster somehow...

Note: See TracTickets for help on using tickets.