Xorshift Generators: Appendix

In his original paper, Marsaglia included (9, 5, 1) as a maximal triplet for xorshift32. That is likely a typo, as the actual maximal triplet is (9, 5, 14).

Using Marsaglia’s triplet (13, 17, 5):

Using Marsaglia’s triplet (13, 7, 17):

• The three shifts cannot all be left shifts, or right shifts (Proposition 4.1)
• If is a maximal triplet, it is also maximal for any of those 8 xorshift variants (Proposition 4.4), because all of their matrix forms are similar
• If is maximal, then is also maximal:

(3)  

This is because their respective transformation matrices are the transposes of each other. Transposing a matrix does not change its characteristic polynomial.

• A triplet of the kind cannot be maximal:

(4)  

This is true because the transition matrix would collapse into:

When a matrix is palindromic (i.e.: ), its characteristic polynomial becomes self-reciprocal (which means that its coefficients read exactly the same way forward and backwards). This makes it non-primitive, hence unable to achieve maximal period.

• For a triplet to be maximal, , , , and have to be setwise coprime:

(5)  

This means they cannot share any prime factor. When is a power of 2, its only prime factor is 2, and this property simply means that at least one among , , and has to be odd (otherwise ).

The following properties are only conjectures, but likely true as were tested all values calculated so far:

• When is even, must be odd:

(6)  

The proposed mechanism for this conjecture is that when and are even, then for every odd . If that happens, then all the odd-degree terms of the characteristic polynomial vanish (because they are linked directly to through Newton’s identities). An even-degree term-only polynomial is a perfect square in characteristic 2, which makes it reducible, disqualifying it from being primitive.

• and cannot be both even at the same time
• A canonical triplet cannot be maximal if . This is likely because the shifts are “pushing out” too many bits, causing too much information to be lost, and collapsing the state.

• Statistical bias torwards small a, small-tomedium b, mid-range c
  • 
  • 
  • 
  • 
• The density of canonical triplets scales like , not
  • Their distribution is not volumetric, and this suggests a lower-dimensional or heavily constrained subset.
  • A possible (unconfirmed) explanation is that the number of monic degree- polynomials is , and the number of monic irreducible degree- polynomials grows asymptotically to . Therefore, a random degree- polynomial is irreducibile with probability . As a consequence, the expected number of maximal triplets among all triplets (which scale with ) would be be approximating .
• Larger correlates with smaller
• Larger correlates with smaller
• and are positively correlated (beyond the canonical constraints )

Multiplicative constant comes from Pierre L’Ecuyer’s research (“Tables of Linear Congruential Generators of different sizes and good lattice structure“)

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

Using Marsaglia’s triplet (11, 8, 19):

uint32_t x, y, z, w;
uint32_t next(void)
{
	uint32_t t = x ^ (x << 11); // A
	x = y;
	y = z;
	z = w;
	return w = w ^ (w >> 19) ^ t ^ (t >> 8); // C, B
}

Unity’s Random class implements a xorshift128 with 4 32-bit words (Unity Documentation: Random).

The C++ source code is not available, but a C# version has been reversed engineered and is available here.

public static uint x = 0, y = 0, z = 0, w = 0;

const uint MT19937 = 1812433253;

public static void InitSeed(int seed)
{
	x = (uint)seed;
	y = (uint)(MT19937 * x + 1);
	z = (uint)(MT19937 * y + 1);
	w = (uint)(MT19937 * z + 1);
}

public static void InitState(uint x, uint y, uint z, uint w)
{
	XORShift128.x = x;
	XORShift128.y = y;
	XORShift128.z = z;
	XORShift128.w = w;
}

public static uint XORShift()
{
	uint t = x ^ (x << 11);
	x = y; y = z; z = w;
	return w = w ^ (w >> 19) ^ t ^ (t >> 8);
}

The seed is initialised using a technique inspired by the Mersenne Twister, using the MT19937 constant 1812433253.

• NVIDIA cuRAND function implements a xorwow160 (CUDA Toolkit Documentation)
• The period is . The official documentation indicates the period is “grater than “
• Uses Marsaglia’s triplet (2, 1, 4), as well as its initial state (which can be found by setting the seed to zero)

uint32_t x = 123456789;
uint32_t y = 362436069;
uint32_t z = 521288629;
uint32_t w = 88675123;
uint32_t v = 886756453;
uint32_t d;

uint32_t next(void)
{
	uint32_t t= x ^ (x >> 2); // A
	x = y;
	y = z;
	z = w;
	w = v;
	v = v ^ (v << 4) ^ t ^ (t << 1); // C, B
	
	return (d+=362437) + v; // Weyl sequence
}

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

• Current implementation used by V8 JavaScript Engine.
• Adopted by Google in 2015 (“There’s Math.random(), and then there’s Math.random()“)
• Uses a different triplet compared to Vigna’s original xorshift128+ implementation
• Implementation code from random-number-generator.h:

uint64_t next (uint64_t* state0, uint64_t* state1) {
	uint64_t s1 = *state0;
	uint64_t s0 = *state1;
	*state0 = s0;
	s1 ^= s1 << 23; // a
	s1 ^= s1 >> 17; // b
	s1 ^= s0;
	s1 ^= s0 >> 26; // c
	*state1 = s1;
	return s0 + s1;
 }

The following techniques rely on left bitwise rotations. Vigna uses in his paper; to avoid confusion with the right bitwise shifts, the letter is used instead.

A left rotation can be implemented using two bitwise shifts:

From which we can extract the following expression:

Hence, the rotation matrix is:

32 bit rotation:

static inline uint32_t rotl(const uint32_t x, int k) {
	return (x << k) | (x >> (32 - k));
}

64 bit rotation:

static inline uint64_t rotl(const uint64_t x, int k) {
	return (x << k) | (x >> (64 - k));
}

• Vigna suggests using (26, 9, 13) for both xoroshiro64* and xoroshiro64**
• The core of Java’s L32X64Mix is based on xoroshiro64 (26, 9, 13)

static uint32_t s[2];

uint32_t next(void) {
	const uint32_t s0 = s[0];
	uint32_t s1 = s[1];
	const uint32_t result = rotl(s0 * 0x9E3779BB, 5) * 5;	// starstar
	const uint32_t result = s0 * 0x9E3779BB;	// star

	s1 ^= s0;
	s[0] = rotl(s0, 26) ^ s1 ^ (s1 << 9); // a, b
	s[1] = rotl(s1, 13); // c

	return result;
}

• Vigna suggests using different triplets for different scrambled variants:
• (24, 16, 37): xoroshiro128**, xoroshiro128+
• (49, 21, 28): xoroshiro128++
• (55, 14, 36): xoroshiro128+ (older version from 2016)

static uint64_t s[2];

uint64_t next(void) {
	const uint64_t s0 = s[0];
	uint64_t s1 = s[1];
	const uint64_t result = s0 + s1;

	s1 ^= s0;
	s[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b
	s[1] = rotl(s1, 37); // c

	return result;
}

• Vigna suggests (25, 27, 36)

static int p;
static uint64_t s[16];

uint64_t next(void) {
	const int q = p;
	const uint64_t s0 = s[p = (p + 1) & 15];
	uint64_t s15 = s[q];
	const uint64_t result = rotl(s0 + s15, 23) + s15;	// plusplus
	const uint64_t result = rotl(s0 * 5, 7) * 9;	// starstar
	const uint64_t result = s0 * 0x9e3779b97f4a7c13;	// star

	s15 ^= s0;
	s[q] = rotl(s0, 25) ^ s15 ^ (s15 << 27);
	s[p] = rotl(s15, 36);

	return result;
}

Vigna reported 42 maximal triplets for xoroshiro2049 (64bit); the correct number is 59.

Vigna reported 25 maximal triplets for xoroshiro4096 (64 bit); the correct number is 37.

It seems that for xoroshiro, when the word size is half the state size, the triplets are symmetrical (i.e.: if is maximal, then is also maximal).

This tested true for xoroshiro32 (16bit word), xoroshiro64 (32bit word), xoroshiro128 (64bit word).

static uint32_t s[4];

uint32_t next(void) {
	const uint32_t result = rotl(s[0] + s[3], 7) + s[0];	// plusplus
	const uint32_t result = rotl(s[1] * 5, 7) * 9;	// starstart
	const uint32_t result = s[0] + s[3];	// plus

	const uint32_t t = s[1] << 9;

	s[2] ^= s[0];
	s[3] ^= s[1];
	s[1] ^= s[2];
	s[0] ^= s[3];

	s[2] ^= t;

	s[3] = rotl(s[3], 11);

	return result;
}

static uint64_t s[4];

uint64_t next(void) {
	const uint64_t result = rotl(s[0] + s[3], 23) + s[0];	// plusplus
	const uint64_t result = rotl(s[1] * 5, 7) * 9;	// starstart
	const uint64_t result = s[0] + s[3];	// plus

	const uint64_t t = s[1] << 17;

	s[2] ^= s[0];
	s[3] ^= s[1];
	s[1] ^= s[2];
	s[0] ^= s[3];

	s[2] ^= t;

	s[3] = rotl(s[3], 45);

	return result;
}

static uint64_t s[8];

uint64_t next(void) {
	const uint64_t result = rotl(s[0] + s[2], 17) + s[2];	// plusplus
	const uint64_t result = rotl(s[1] * 5, 7) * 9;	// startstar
	const uint64_t result = s[0] + s[2];	// plus

	const uint64_t t = s[1] << 11;

	s[2] ^= s[0];
	s[5] ^= s[1];
	s[1] ^= s[2];
	s[7] ^= s[3];
	s[3] ^= s[4];
	s[4] ^= s[5];
	s[0] ^= s[6];
	s[6] ^= s[7];

	s[6] ^= t;

	s[7] = rotl(s[7], 21);

	return result;
}