COIN-OR::LEMON - Graph Library

Opened 16 years ago

Last modified 13 years ago

#105 new enhancement

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

Reported by: Alpar Juttner Owned by: Alpar Juttner
Priority: major Milestone:
Component: core Version: hg main
Keywords: random Cc: Balazs Dezso
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 Szabó Mátyás 14 years ago.
Ziggurat.mw (113.8 KB) - added by Szabó Mátyás 13 years ago.
Maple sheet I used to generate the constants for random.cc
3ae5695818c3-b33813f266e3.patch (71.0 KB) - added by Alpar Juttner 13 years ago.
Virtual Campus Supercomputing Center.png (14.0 KB) - added by Smartmil8 9 years ago.
about.me purevolume.com
developedBy.jpg (4.0 KB) - added by Smartmil8 9 years ago.
wardsauto.com winsupersite.com

Download all attachments as: .zip

Change History (17)

comment:1 Changed 15 years ago by Alpar Juttner

Milestone: LEMON 1.1 release

Changed 14 years ago by Szabó Mátyás

Attachment: 2fefbb0c8a28.patch added

comment:2 Changed 14 years ago by 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).

comment:3 Changed 13 years ago by Alpar Juttner

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

Changed 13 years ago by Szabó Mátyás

Attachment: Ziggurat.mw added

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

comment:4 Changed 13 years ago by Szabó Mátyás

I hope it's understandable. :)

Changed 13 years ago by Alpar Juttner

comment:5 Changed 13 years ago by 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.

comment:6 Changed 13 years ago by Alpar Juttner

Cc: Balazs Dezso added
Priority: minormajor

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 13 years ago by 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 :)

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

comment:9 in reply to:  8 Changed 13 years ago by 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.)

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

comment:11 in reply to:  10 Changed 13 years ago by 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 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 13 years ago by 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 231. So we can multiply our consts by 231 (and by 2(-31)), and use them instead. Which makes it faster somehow...

Changed 9 years ago by Smartmil8

Attachment: developedBy.jpg added
Note: See TracTickets for help on using tickets.