static unsigned long Q[4096], c=123;
unsigned long CMWC(void)
{
unsigned long long t, a=18782LL;
static unsignd long i=4095;
unsigned long x, m=0xFFFFFFFE;
i = (i + 1) & 4095;
t = a * Q + c;
c = (t >> 32);
x = t + c;
if(x < c)
{
++x;
++c;
}
return (Q = m - x);
}
Number Generator "CMWC 4096" and Implementation?
Author
In another thread (C++ true random numbers), somebody mentioned that the generator CMWC 4096 by George Marsaglia was better in several ways (speed, period, simplicity, output quality) than the more commonly used Mersenne Twister. These details spurred me into investigating this generator, and I need some assistance doing so.
A google quickly produced (not-so-great) C source code created by the creator of the CMWC algorithm (below). You can see a paper by the same person (taken from wikipedia) titled Random number generators. An explanation of the algorithm is also available on wikipedia: Complementary-multiply-with-carry RNGs.
I think I understand the premise of the algorithm after reading the paper and several other comments by the author, but there are a couple of things in the code that don't seem to match the explanations.
(1) To me, it looks like the carry is being added twice (once to 't' and again when the lower bits are extracted to 'x'). I would believe this is simply a basic mistake because it doesn't match my understanding of the math, but every reproduction of the algorithm seems to contain the apparent double addition of 'c'.
(2) I also do not understand the purpose of the "if" statement and it's body - I don't see anything corresponding to it in his very simple explanation of the algorithm in the paper I linked above. I imagine that if (1) is not a mistake that this is related to it (testing for overflow), but I'm not aware of the reason for (1) so this also remains a mystery to me.
Author
I realized something before I posted the above, but I apparently forgot while making the post. With respect to (1), the carry isn't added twice, but rather both Carry_n-1 and Carry_n are added. I still don't have any idea why, but it is an important distinction that two different carry values are added.
This does not answer your question, but one of the first google results presents an implementation given by the author, and he includes a
if ((x+1)==0) { c++; x=0; }
between the if (x < c) and the return.
I just figured it might be relevant since it includes the carry variable as well.
I don't understand why he does that though... doesn't that knock out the only chance that the function will return 0xffffffff?
-Scott
if ((x+1)==0) { c++; x=0; }
between the if (x < c) and the return.
I just figured it might be relevant since it includes the carry variable as well.
I don't understand why he does that though... doesn't that knock out the only chance that the function will return 0xffffffff?
-Scott
I agree. That code does not calculate the same function as documented on the wikipedia page. His cn is calculated as (a * Q + cn-1) mod b div b, and his x is calculated as (a * Q + cn-1) mod b + cn mod b then biased by 1 if x < cn (where b == 2^(2^sizeof unsigned long) and sizeof unsigned long == 4).
I can't for the life of me come up with any kind of reasonable algebraic transformation from the wikipedia equation and the implementation. Perhaps there are some cool manipulative tricks I've missed (such as Schrage's algorithm, only different).
I would also argue that a 16 kilobyte state vector per RNG is pretty expensive compared with the MT (or even other lag generators).
I can't for the life of me come up with any kind of reasonable algebraic transformation from the wikipedia equation and the implementation. Perhaps there are some cool manipulative tricks I've missed (such as Schrage's algorithm, only different).
I would also argue that a 16 kilobyte state vector per RNG is pretty expensive compared with the MT (or even other lag generators).
I'll give it a try:
static unsigned long Q[4096], c=123; // Q holds the stateunsigned long CMWC(void){ unsigned long long t, a=18782LL; static unsignd long i=4095; unsigned long x, m=0xFFFFFFFE; i = (i + 1) & 4095; // Advance through the state array, circularly t = a * Q + c; // Notice that the result is computed as a 64-bit integer. // This is the core of a traditional LGC, c = (t >> 32); // except we keep changing the constant we add. x = t + c; // add the lower 32 bits with the higher 32 bits to avoid poor quality bits at the extremes. if(x < c) // About 50% of the time (?? I don't know why this helps). { ++x; ++c; } return (Q = m - x); // I don't know why this is better than simply `(Q = x)'.}
Author
Quote:Original post by BregmaIn some applications, such a large state is required. For example, consider a card game where each player can have a deck up to size S and there can be N players per game, you need to be able to generate a number from 0 to S!^N which grows very fast. If the game allowed up to 1000 cards per deck and 6 players per game, the state size required is 1000! ~= 2^8520 for one deck and (2^8520)^6 = 2^51120 for the initial game state, which is beyond MT's capability.
[...]I would also argue that a 16 kilobyte state vector per RNG is pretty expensive compared with the MT (or even other lag generators).
This kind of situation could easily happen in a customizable card game that doesn't have a maximum deck size in the rules, such as Magic the Gathering. In fact, to allow some tournament formats, you'd actually need two levels of RNGs (probably one 'seed generator' per tournament and then one 'game generator' per game) since a single one does possible provide the required state size.
Quote:I'll give it a try:It's not 50% of the time: Since C is added to the value that determines X, X can only be lower than C if the value has overflowed. Overflow can only happen rather rarely since the range of C (0..18782 due to the ++C) is rather low and the range of the low 32 bits of T (0..2^32-1) is rather high. 99.99956272% of the time ((2^32-18782)/(2^32-1)), the lower bits of T will be small enough that no value of C will cause overflow.
[...]if(x < c) // About 50% of the time (?? I don't know why this helps).[...]
[Edited by - Extrarius on October 22, 2008 11:22:48 AM]
Quote:Original post by ExtrariusQuote:I'll give it a try:It's not 50% of the time: Since C is added to the value that determines X, X can only be lower than C if the value has overflowed. Overflow can only happen rather rarely since the range of C (0..18782 due to the ++C) is rather low and the range of the low 32 bits of T (0..2^32-1) is rather high. 99.99956272% of the time ((2^32-18782)/(2^32-1)), the lower bits of T will be small enough that no value of C will cause overflow.
[...]if(x < c) // About 50% of the time (?? I don't know why this helps).[...]
Of course you are right. When I first read the algorithm, I thought that Q was a 64-bit integer.
Anyway, it's not clear to me what that `if(x < c)' part does, and the same goes for the subtracting-from-0xfffffffe part.
Quote:Original post by alvaro
and the same goes for the subtracting-from-0xfffffffe part.
The same as subtracting 4 bit number.
14 - 0 = 14
14 - 1 = 13
14 - 15 = -1 mod 0xF = 15
Basically inverting the number, and shifting it by a little.
Complement Multiply With Carry
I am not capable of translating the theory to an algorithm either, but I want to point out that the calculation of Mod (2^n-1) needs to use some tricks.
For example, if modulus = 2^n - 1 then there is an efficiency trick:
result = x Mod modulus
is equivilent to
result = (x And modulus) + (x >> n)
If (result >= modulus) Then result -= modulus
(^^ thats a bitwise And)
The later is generally much more efficient (assuming good branch prediction)
Now, note that he is using the constant 0xfffffffe, which is not his modulus (which is 0xffffffff) but in fact 1 less than his modulus .. this could explain what might be considered a spurious increment
note that (x And modulus) in his code is implicit with a conversion from 64-bit to 32-bit
I am not capable of translating the theory to an algorithm either, but I want to point out that the calculation of Mod (2^n-1) needs to use some tricks.
For example, if modulus = 2^n - 1 then there is an efficiency trick:
result = x Mod modulus
is equivilent to
result = (x And modulus) + (x >> n)
If (result >= modulus) Then result -= modulus
(^^ thats a bitwise And)
The later is generally much more efficient (assuming good branch prediction)
Now, note that he is using the constant 0xfffffffe, which is not his modulus (which is 0xffffffff) but in fact 1 less than his modulus .. this could explain what might be considered a spurious increment
note that (x And modulus) in his code is implicit with a conversion from 64-bit to 32-bit
Quote:Original post by RagharQuote:Original post by alvaro
and the same goes for the subtracting-from-0xfffffffe part.
The same as subtracting 4 bit number.
14 - 0 = 14
14 - 1 = 13
14 - 15 = -1 mod 0xF = 15
Basically inverting the number, and shifting it by a little.
When I say I don't know what the code does I mean that I don't see how it's going to significantly affect the quality of the pseudo-random numbers generated.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
