Jump to content
  • Advertisement
Sign in to follow this  
Extrarius

Unity Number Generator "CMWC 4096" and Implementation?

This topic is 3560 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

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



Share this post


Link to post
Share on other sites
Advertisement
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.

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
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).

Share this post


Link to post
Share on other sites
I'll give it a try:
static unsigned long Q[4096], c=123; // Q holds the state
unsigned 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)'.
}


Share this post


Link to post
Share on other sites
Quote:
Original post by Bregma
[...]I would also argue that a 16 kilobyte state vector per RNG is pretty expensive compared with the MT (or even other lag generators).
In 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.

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:
[...]if(x < c) // About 50% of the time (?? I don't know why this helps).[...]
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.

[Edited by - Extrarius on October 22, 2008 11:22:48 AM]

Share this post


Link to post
Share on other sites
Quote:
Original post by Extrarius
Quote:
I'll give it a try:
[...]if(x < c) // About 50% of the time (?? I don't know why this helps).[...]
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.


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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
Quote:
Original post by Raghar
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.


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.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!