Xorshift Generators

This article is a deep dive into xorshift generators, one one the most adopted family of algorithms to generate pseudorandom numbers. This is a companion to a video documentary which will be released in the following weeks. If you want to learn more about xorshift generators, including tables of maximal triplets, you can consult the Xorshift: Appendix. Introduction At this very moment, billions of devices are running these three lines of code. uint64_t x; uint64_t y; uint64_t next() { uint64_t t = x; const uint64_t s = y; x = s; t ^= t << 23; // a t ^= t >> 17; // b t ^= s ^ (s >> 26); // c y = t; return t + s; } They hide inside virtually every browser, and the servers they talk to.They control the events that happen while you’re playing. From the worlds you explore, to the cards in your hand. This algorithm—and the others like it—are the metaphorical dice that your computer rolls every time it makes a decision. But …what does it do, and how does it work? I’m Alan Zucconi, and in this documentary we’ll travel to the very source of randomness, by dissecting the most popular family of pseudorandom algorithms: xorshift generators. And we’ll do it by asking a very simple question: what makes 23, 17 and 26 special? Part 1: Xorshift Generators Pseudorandom numbers Computers are terrible at being random. In fact, every “random” number you have ever generated was almost certainly completely deterministic. But it turns out that with the right algorithm, deterministic sequences can look random enough to power games, simulations, procedural content, cryptography—well, okay, maybe not cryptography …but we’ll get back to that later. Most of these pseudorandom techniques work on the same principle: there’s a function that we apply to the previous number in the sequence, to get the next one. The first one—the seed—is given, and the rest of the sequence follows from it. That’s how games like No Man’s Sky and Minecraft are able to “compress” entire galaxies and worlds in a number. That’s the seed that kickstarts the sequence of pseudorandom numbers, and with the same seed, you get the same sequence. If you are curious to find out more about this, I made an entire documentary about the world generation of Minecraft. Marsaglia’s Xorshift In 2003, mathematician George Marsaglia published a short paper (“Xorshift RNGs“) introducing a family of generators that fall into this category, based almost on two operations: XOR and shifts. Unsurprisingly, they’re called xorshift generators. This one here, for example, generates sequences of 32-bit numbers: uint32_t x; uint32_t xorshift32(void) { x ^= x << 13; x ^= x >> 17; x ^= x << 5; return x; } 📝 Variable types used by Marsaglia In his original paper, Marsaglia’s code used unsigned long for the variables. In C/C++, the size of types depends on the platform. In his case, unsigned long is 32 bits, and unsigned long long is 64 bits. To avoid confusion, in this article we used uint32_t and uint64_t. If you are implementing this in C#, you can use uint and ulong, respectively. ✏ Correction The original paper published by Marsaglia contained a typo. The xorshift32 code on page 4: unsigned long xor(){ static unsigned long y=2463534242; yˆ=(y<<13); y=(y>>17); return (yˆ=(y<<5)); } Should have been instead this: unsigned long xor(){ static unsigned long y=2463534242; yˆ=(y<<13); y^=(y>>17); return (yˆ=(y<<5)); } 📝 Variants There are 8 possible variants of this code: ACodeXx^= x << a; x ^= x >> b; x ^= x << c;x^= x >> a; x ^= x << b; x ^= x >> c;x^= x << c; x ^= x >> b; x ^= x << a;x^= x >> c; x ^= x << b; x ^= x >> a;x^= x << a; x ^= x << c; x ^= x >> b;x^= x >> a; x ^= x >> c; x ^= x << b;x^= x >> b; x ^= x << a; x ^= x << c;x^= x << b; x ^= x >> a; x ^= x << c;xorshift32 and xorshift64 variants, as classified by Vigna (A column) and Panneton and L’Ecuyer (X column). Table from this paper. If a triplet is maximal for one of those, it is also maximal for the other seven (although the overall sequence might exhibit different statistical properties). Let’s see how it works, line by line! It takes the previous number in the sequence (x), makes a copy, and shifts its bits 13 times to the left. It drops the bits that overflow, and fills the rest with zeros. left shift x << 1 right shift x >> 1 Now: every time there’s a “1”, it flips the corresponding original bit: that’s a bitwise XOR. Next, it does another xorshift: copy, shift, XOR. But this time, shifting by 17 bits to the right. And once more, by 5 bits to the left. Diagram of xorshift32 No multiplications. No divisions. No lookup tables. Just a few bitwise instructions. If we paint a texture using these values as colours …well, it looks pretty random to me! xorshift32 (13, 17, 5) But why did Marsaglia himself recommend shifting by 13, 17, and 5 bits? Why not 42, 67, and 1729? And why exactly three xorshift operations, and not two, or four? It’s questions like these that you should always ask when you see magic numbers in someone else’s code. Like 0x5F3759DF from Quake’s Fast Inverse Square Root function, or the golden ratio constant 0x9e3779b97f4a7c15 found in a lot of hashing algorithms. 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. 📝 Period length When a triplet is maximal, the sequence of pseudorandom numbers it generates cycles through all possible non-zero values. This means that changing the seed effectively “shifts” where the sequence starts from. When the triplet is not maximal, the period is partitioned across smaller cycles. Different seeds will land in one of the many cycles of the generator, but none of them will be maximal. You can read more about the underlying mathematics in this paper titled “The cycle structure of a linear transformation over a finite field“. 📝 Canonical triplets For “monolithic” xorshift generators only, Marsaglia introduced the concept of canonical triplets. A canonical triplet is a triplet where . This is relevant, because: If is maximal, then is also maximal. This can help reduce the search space. A triplet of the form cannot be maximal. This means that there are twice as many maximal triplets as there are maximal canonical triplets. All of these results are only true for “classical”, monolithic xorshifts; for other techniques (such as Marsaglia’s companion block matrices, xoroshiro, xoshiro, …) this is not necessarily true, and the introduction of canonical triplets is irrelevant. 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! shifts \ bits32 bits64 bits2 shifts<< a>> bno solutionsxorshift64(7, 9)2 solutions3 shifts<< a>> b<< cxorshift32(13, 17, 5)162 solutionsxorshift64(13, 7, 17)550 solutionsTable of popular xorshift algorithms using 2 and 3shifts ✏ Correction In his original paper, Marsaglia indicated there were no maximal tuples for 64 bits. That is not the case, as both (7,9) and (9,7) are maximal. 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 . ❓ Primitive lines You can clearly see in the scatter plot above a “line”, emerging from a sequence of aligned maximal triplets. All points on that line have different matrices which share the same primitive polynomial. Such “primitive lines” appears to lie on the “edges” and “faces” of the three-dimensional structure that emerges from the maximal triplet distribution. They becomes really abundant when the period is a Mersenne prime. 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… ❓ The Hull-Dobell Theorem Finding the maximal period is not a problem unique to xorshift generators. Linear Congruential Generator (LCG) are another type of well-known pseudo-random generators, governed by the following recurrence: (1) The Hull-Dobell Theorem provides the exact mathematical conditions required for a LCG to achieve its maximum possible period: The modulus and the increment must be relatively prime (their greatest common divisor is 1) The multiplier must be divisible by every prime factor of the modulus If the modulus is divisible by 4, then must also be divisible by 4 When those conditions are met, the LCG has maximal period. For xorshift generators, things are unfortunately more complex… Part