Last time, we saw how to use rand() and why rand() is unacceptable. We also took a look at the Mersenne Twister, which is definitely a huge improvement over rand(), but we were put off by its memory footprint of 2.5K. So, how can we do better than Mersenne Twister?

Enter George Marsaglia, in a posting to Sci.Stat.Math on January 12th, 1999. (Either that, or December 1st, 1999, depending on how Bob Wheeler writes his dates.) In that post, he puts down 17 lines of code that implement 8 different RNGs. That sounds really promising, at least as a starting point. Let’s pick one of those and run with it and see where it takes us.

The RNG that we’re going to be using as our basis is Multiply With Carry (MWC). Let’s isolate Marsaglia’s code for just the MWC RNG:

#define UL unsigned long
#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
#define MWC   (znew+wnew)
static UL z=362436069, w=521288629;
void settable(UL i1, UL i2)
 { z=i1; w=i2; }

Hrm. The first obvious problem is that this RNG uses a static context. Fortunately, it’s small and simple enough that we can fix that on our own by taking z and w and putting them into a struct or class. But before we go down that rabbit hole, let’s examine what this code is actually doing:

We’ll start with the MWC macro. MWC adds the newly-generated Z value to the newly-generated W value. Fair enough. What do those do? The first thing you’ll notice is that there’s a pattern to the Z and W updates, so we’ll take apart the algorithm for the Z/W update:

  • Split off the 32-bit value into two 16-bit values, which we’ll call L for low and H for high.
  • Multiply L by the constant C (36969 for Z and 18000 for W).
  • Add the result to H. This is your new Z or W value.
  • Return the bottom 16 bits of the value. For Z, shift them over into the high 16 bits of the return value.

From this, we can determine that Z and W are 32-bit values, each of which represents a 16-bit RNG. The MWC RNG generates its 32-bit value by combining the two 16-bit RNGs into a single value in the high and low words of the value. Now that we understand that, what does it mean? It means that we have to pick *good* values for our initial Z and W, because if we don’t, then we really only have 16 bits worth of randomness.

It kind of sucks to have to pick two seeds for an RNG. Can we improve on that? How can we take Marsaglia’s original MWC RNG and improve it by just requiring one seed? Well, let’s look at what we determined about the Z and W values: each one is a 32-bit value that represents a 16-bit portion of the 32-bit RNG. Maybe we can extend that concept by selecting a 64-bit value which represents a 32-bit RNG. We’ll start by just updating Marsaglia’s code, then examine the process of making it practical for every-day use:

#define ULL unsigned long long
#define vnew  (v=36969*(v&0xFFFFFFFF)+(v>>32))
#define MWC   (vnew&0xFFFFFFFF)
static ULL v=0x159A55E51F123BB5; // 64-bit value compiled from original z and w as z<<32 + w
void settable(ULL i1)
 { v=i1; }

Beautiful! It’s tiny, it’s fast, it’s deterministic, it’s got decently long cycles…the only problem is the static context. Fortunately, that one is easy to fix. Let’s start by making a structure to hold the V value:

struct RAND
{
	unsigned long long v;
};

Now a function to initialize it:

#define DEFAULT_V_VALUE 0x159A55E51F123BB5
inline void RandInit(RAND& tRand, unsigned long long seed)
{
	if(seed == 0)
		seed = DEFAULT_V_VALUE;

	tRand.v = seed;
}

And finally, a function to actually generate a random number:

inline unsigned int RandGetNum(RAND& tRand)
{
	tRand.v = 36969 * (tRand.v & 0xFFFFFFFF) + (tRand.v >> 32);
	return tRand.v & 0xFFFFFFFF;
}

So there you have it! A small, fast, deterministic, long-cycled, non-static context RNG.

Next time: Exploring Compiler Optimizations