These constants might be arbitrary, but they’re never random. Because they almost always reveal something deeper about the problem, and the way in which it was solved.

So… what happens if we try a different triplet!?

These sequences don’t look that random anymore! And that’s because the period of most triplets—which is, how many unique numbers they generate—is far too small to be used for procedural generation!

What makes 13, 17, and 5 special is that it’s a maximal triplet! Out of all the possible ones—that’s —there are only 162 that generates all 32-bit numbers. Well, all except zero, which xorshifts can never generate anyway.

Over four billion states () from three lines of code.

If you’re thinking this works because these (13, 17, and 5) are all prime numbers… well, you’d be wrong! In fact, (2, 15, 25) and (16, 21, 9) work… just as fine!

And if you test all combinations, you’ll also find out that we need 3 xorshifts, because you can’t generate all 32-bit numbers with just two of them!

When visualised in three dimensions, the maximal triplets arrange on a complex structure. Below, you can see the arrangement of the canonical maximal triplets for .

As gets larger, the structure appears to get denser, not richer. This indicates that the lattice structure of maximal triplets is consistent across different values of .

Finding triplets

But how do we find those maximal triplets? The most naïve way is to count how many unique numbers each one generates. With 32 bits we have almost 30.000 triplets to check (), and each one generates up to almost 4.3 billion unique numbers (). Which is… a lot, but nothing we can’t test in a few hours!

So what about 64-bit numbers? You may think that doubling the number of bits, also doubles the search time. But that would be …incorrect! It is every new bit that doubles the range of numbers we can represent. So going from 32 to 64 bits actually doubles the range—and the search time—32 times! That’s roughly …49 million years. Give or take.

So how could Marsaglia do this—in minutes—back in 2003? Well, he did it with two very clever tricks…

Part 2: Matrix Approach

Underneath their code, xorshift generators are secretly just matrix multiplications:

(2)  

And this matters, because Linear Algebra can predict the length of a sequence, without calculating all of its elements.

Let’s treat the bits of our seed as the components of a vector. Every instruction in our code translates to a linear transformation between vectors and matrices. That’s our first trick, and it’s not a coincidence: the algorithm is literally designed like this to make the search for maximal triplets feasible.

Bitshifts

You probably already know that when a vector is multiplied by the identity matrix, it remains unchanged.

(3)  

This happens because matrix multiplication is defined in such a way that the positions of the 1s in each column define which bits are linearly combined to produce the output vector.

Shifting the diagonal down shifts the output bits to the left. Let’s call this matrix !

(4)  

If we want to shift by two bits (x << 2), we multiply the result by once more!

And because multiplication is associative, it turns out that to shift a vector by bits (x << a), we multiply it by raised to the -th power.

(5)  

And the exact same reasoning holds for a right shift (x >> a), which is equivalent to :

(6)  

Xor

The bitwise XOR (between two variables) happens to behave exactly like the component-wise addition (between two vectors).

And since all those vectors contain either zero or one, let’s say that 1+1=0. You can think about this as “addition without carry over”, or “addition followed by modulo 2”. But—for the hardcore fans of Field Theory—that’s simply a finite field with two elements.

Xorshift

We now have everything we need to express one single xorshift operation (x ^ (x << a)) in its matrix form: .

And if we want to perform three of them in series, we just multiply the original number by three separate matrices!
That’s our final matrix , and it looks like this:

(7)  

It is important to remember that those are binary matrices: every cell can only contain zero or one. This means that we can visualise them as bitmaps:

And when we do that, the matrix looks like this:

Part 3: Order Test

Now that the xorshift is framed in terms of vectors and matrices, we can take advantage of Linear Algebra. The period of a generator—which is, how many unique numbers it produces—corresponds to the multiplicative order of its matrix.

That is defined as the smallest value for which :

(8)  

Hence, a triplet is maximal if and only if that multiplicative order is equal to all the possible non-zero numbers that can be represented with bits:

(9)  

Let’s see what all of this means in practice.

Each time we multiply by , we progress in the sequence.

(10)  

which means that:

(11)  

Matrix multiplication is associative, so we can group all the s together. The -th element of the sequence is found by multiplying the first one by , skipping all the intermediate ones!

(12)  

After a certain number of iterations—let’s call it for period—we get back to the starting element, and the sequence repeats:

(13)  

Since is guaranteed to be non-null (as xorshift generators can never produce 0), we can simplify the s on both sides, and find out that if a sequence repeats after iterations, then is the identity matrix! Let’s call this the identity test.

Maximal Order test

Now, if a number p passes the identity test, it doesn’t necessarily mean that’s the actual period of our sequence! If numbers repeat after 30 iterations; they’ll also repeat after 60, after 90, after 120 …because it’s like looping through the sequence two, three, four times…

In a nutshell, if passes the identity test, so do all of its integer multiples (with ):

(14)  

This also means that if 30 passes the identity test, perhaps that’s because the actual period is 15 and we just looped twice! Or perhaps it was 10, and we looped 3 times. Or it was 6, and we looped 5 times!

You get where I’m going with this: for 30 to be the smallest period, none of its divisors should pass the identity test—including 5, 3, and 2.

In fact, (\eq{Tn_eq}) also implies that if a number doesn’t pass the identity test, then none of its proper divisors will:

(15)  

But we don’t have to test them all! 15 is 5 times 3: if 15 doesn’t pass the identity test, neither 5 nor 3 do! If 10 doesn’t pass the identity test, neither 5 nor 2 do. And the same for 6!

So what’s the minimum number of divisors to test?

Well, the prime factors of 30 are 2, 3, and 5; and its divisors are combinations of them.
Testing larger numbers first helps prune more divisors out, so we can remove one prime factor at a time from 30, and those are the minimum number of tests that covers all other cases.

Going back to our original 64bit problem: for a triplet to have period , itself must pass the identity test, while that same test must fail when we divide by each of its prime factors.

That reduces the number of iterations from 127–all possible divisors including itself–to just 8!

That’s the second trick! And that’s how we find that with 64 bits, we have 550 maximal triplets. With Marsaglia choosing 13, 7 and 17.

Part 4: Companion Matrices

We’ve done 32 bits… we’ve done 64… what about 128?

Our algorithm finds 1072 full-period triplets, but none of them can be found in the selected list presented by Marsaglia: (11, 8, 19), (5, 14, 1), (15, 4, 21), (23, 24, 3), and (5, 12, 29).

What’s going on? Well, modern architectures natively support bitwise operations with 32 and 64 bits. But that’s not typically the case for 128, making a “true”, “monolithic” xorshift128 neither efficient to run, nor elegant to write.

uint32_t x;
uint32_t y;
uint32_t z;
uint32_t w;

uint32_t next()
{
	uint32_t tx, ty, tz, tw;

	// s ^= s << a
	tx = (x << a) | (y >> (32-a));
	ty = (y << a) | (z >> (32-a));
	tz = (z << a) | (w >> (32-a));
	tw = (w << a);

	x ^= tx;
	y ^= ty;
	z ^= tz;
	w ^= tw;

	// s ^= s >> b
	tx = (x >> b);
	ty = (y >> b) | (x << (32-b));
	tz = (z >> b) | (y << (32-b));
	tw = (w >> b) | (z << (32-b));

	x ^= tx;
	y ^= ty;
	z ^= tz;
	w ^= tw;

	// s ^= s << c   
	tx = (x << c) | (y >> (32-c));
	ty = (y << c) | (z >> (32-c));
	tz = (z << c) | (w >> (32-c));
	tw = (w << c);

	x ^= tx;
	y ^= ty;
	z ^= tz;
	w ^= tw;

	return w;
}

To get around this, Marsaglia proposed a similar but different algorithm: one that still uses bitwise operations, but stores the 128-bit state in four 32-bit integers: x, y, z, and w:

uint32_t x;
uint32_t y;
uint32_t z;
uint32_t w;

uint32_t next()
{
	uint32_t t;
	t = x ^ (x << 11); // a
	t = t ^ (t >>  8); // b
	x = y;
	y = z;
	z = w;
	w = w ^ (w >> 19); // c
	w = w ^ t;

	return w;
}

unsigned long xor128(){
	static unsigned long x=123456789,y=362436069,z=521288629,w=88675123;
	unsigned long t;
	t=(xˆ(x<<11));x=y;y=z;z=w; return( w=(wˆ(w>>19))ˆ(tˆ(t>>8)) );
}

Let’s see how this version works!

Lines 10-12 shift the lower 96 bits to the right, while the rest reshuffles x and w into the lower 32 bits.

The auxiliary variable t stores the result of a left xorshift on x, followed by a right xorshift.

w goes through a right xorshift, and the result is xor-ed with t.

We can work backwards to reconstruct the matrix that maps the input state to the output state and a linear combination of and .

And that’s what it looks like, both in block matrix form …and in its full glory.

(20)  

This looks dangerously similar to a Frobenius companion matrix—studied in the context of primitive polynomials for their connection to maximal sequences—hence making it a “natural choice” to restrict the search.

Part 5: Characteristic Polynomials

Later papers explored that connection further, ditching Marsaglia’s naïve order test, in favour of more targeted approaches. After all, we don’t need to find all maximal triplets: we just need to find one that’s good. And primitive polynomials have played a key part in that.

Every square matrix is paired with a characteristic polynomial, which is defined as:

(21)  

The symbol is often used for the variable (instead of the more common ), especially when eigenvalues and eigenvectors are involved. Regardless, the result is a polynomial of degree , which is used to study certain properties of the matrix.

As it turns out, there’s a special connection between the maximality of a triplet, and its characteristic polynomial: its primitivity.

A polynomial is primitive over if any of its roots is a primitive element of . In a nutshell, it means that:

(22)  

It means that integer powers of a root are cycling through all non-zero elements. Which is exactly what the multiplicative order or the matrix captured.

Krylov Space

Testing for the primitivity of the characteristic polynomial is, generally speaking, much more efficient than testing the multiplicative order of the matrix. Matrix calculations are, in fact, expensive: the multiplication of a dense matrix by a vector has a computational complexity of , while multiplying two dense matrices is .

A “naïve” test for primitivity would first build the matrix , then extract the characteristic polynomial by solving . All of those steps are very expensive, and would likely counteract any actual speedup when compared to the traditional order test seen before.

What makes this approach feasible, is that we can reconstruct the characteristic polynomial without building its associated matrix . We do that by constructing the so-called Krylov space first, a part of which is then used by the Berlekamp-Massey algorithm to reconstruct the characteristic polynomial.

A Krylov space is built by multiplying a starting vector by a matrix over and over:

(23)  

Normally, this would be expensive as it requires not only the construction of , but also performing a vector multiplication for every element of the space.

For xorshift generators, that is not necessary: the Krylov sequence is nothing more than the sequence of states generated by the pseudorandom algorithm, which is exceptionally fast to compute!

Berlekamp–Massey Algorithm

The sequence is then passed to the Berlekamp–Massey algorithm, an algorithm which reconstructs minimal polynomials from linear recurrent sequences in arbitrary fields.

The minimal polynomial is not exactly what we’re after, but it’s close! It’s guaranteed to be a divisor of the characteristic polynomial.

From this, we find ourselves in two scenarios:

This means that the BM algorithm either returns the characteristic polynomial, or one of its proper factors (revealing that the original triplet was not maximal). What it cannot do, on its own, it so tell us if a triplet is maximal. We still don’t know if the returned characteristic polynomial is irreducible, let alone primitive.

But once extracted, the characteristic polynomial can be tested for both conditions, using the same order test previously devised for matrices (but applied to polynomials over ):

Frobenius Companion Matrices

We have learnt in the previous section that a triplet is maximal if its associated characteristic polynomial is primitive. This connection goes both ways: every primitive polynomial represents a linear recurrence from which a generator can be extracted.

Rather than finding maximal triplets by brute force testing their maximality, it is sometimes more efficient to find a polynomial that is known to be primitive, and rebuild its associated matrix.

The simplest way to do that is through Frobenius companion matrices, which we have briefly mentioned before:

(24)  

whose characteristic polynomial is expressed in the following form:

(25)  

The beauty of a companion matrix is that it contains the coefficients of its characteristic polynomial. And if that is known to be primitive, its companion matrix will be maximal.

This is one of the ideas behind Marsaglia’s companion block matrices, which were indeed designed with a similar block structure in mind.

Part 6: Testing

Just because a triplet is maximal, it doesn’t mean it’s also… good.

Different full-period sequences have different statistical properties, and not all of them make good random number generators.

When a sequence looks predictable, it produces linear artefacts which can manifest pretty strongly in procedural terrains and textures.

uint64_t x;
uint64_t next(void)
{
	x ^= x >> 12; // a
	x ^= x << 25; // b
	x ^= x >> 27; // c
	return x * UINT64_C(2685821657736338717);
}

uint64_t x;
uint64_t y;
            
uint64_t next()
{
	uint64_t result = x + y;
	uint64_t t;
	t = x ^ (x << 23); // a
	t = t ^ (t >> 18); // b
	x = y;
	y = y ^ (y >>  5); // c
	y = y ^ t;
	return result;
}

uint64_t s[2];
uint64_t next(void)
{
	uint64_t s1 = s[0];
	const uint64_t s0 = s[1];
	const uint64_t result = s0 + s1;
	s[0] = s0;
	s1 ^= s1 << 23; // a
	s[1] = s1 ^ s0 ^ (s1 >> 18) ^ (s0 >> 5); // b, c
	return result;
}