Mersenne Twister

Every time you call random.random() in Python, a generator called MT19937 runs, shuffles some bits around, and hands you back a number. MT19937 is an implementation of the Mersenne Twister, a pseudorandom number generator built in 1997 by Makoto Matsumoto and Takuji Nishimura.

It’s named after Mersenne primes, primes of the form $2^p - 1$, because the generator’s cycle length before it repeats is exactly $2^{19937}-1$, one such prime. It’s the engine behind Python’s random module, and most other languages use some version of it too.

Each of those numbers depends on a starting value called a seed. Most of the time nobody sets one directly though, Python just picks one for you automatically, by asking the operating system for some bytes of its own cryptographically secure randomness.

Python itself doesn’t generate this randomness, it just delegates to whatever secure source the OS provides: Linux uses a ChaCha20-based generator built into the kernel, macOS uses its own kernel CSPRNG, and Windows uses Microsoft’s CryptGenRandom.

That’s a whole topic on its own, so for this post we’ll skip it and just assume you called random.seed(42) yourself.

Setup

The Mersenne Twister can’t start from nothing. It needs a list of starting numbers to fill its internal memory before it can do anything.

Python’s internal memory holds 624 numbers, each 32 bits wide. Your seed needs to become a list of numbers that Python can use to fill that memory. This list is called the key.

For a small seed like 42, this is simple, the key is just that one number, wrapped in a list:

For a much bigger seed, Python would split it into multiple 32-bit chunks instead, but the idea is the same. Either way, we now have a short list of numbers ready to feed into the actual algorithm. Getting this list wasn’t the Mersenne Twister itself, it’s just Python figuring out what to feed it.

Seeding

Before your actual key gets used, Python fills its 624-slot memory with a fixed, repeatable pattern. This step always starts from the exact same number, 19650218, no matter what your seed is. Think of it as priming the pump before pouring in your specific ingredients.

The rule for filling each slot from the one before it is:

Here $f = 1812433253$ is just a fixed constant, and $\oplus$ means XOR (flip bits where they differ).

Worked example. The first couple of slots, filled in by hand:

This repeats for all 624 slots. Notice your seed hasn’t even been used yet, this part is always identical no matter what.

Mix seed

Now Python walks through the memory a second time, this time actually folding your key into each slot:

Since our key only has one number in it (42), this step reuses that same number over and over, once for each of the 624 slots.

Worked example. The very first mixing step, updating slot 1:

This repeats across all 624 slots, followed by one more short cleanup pass with a different multiplier. After all of it finishes, three of the resulting numbers, the ones we’ll need in a moment, look like this:

Twist

This next step is where the algorithm gets its name, and where its famously long, never-repeating cycle comes from. To generate a new number for a given slot, three of the current numbers get mixed together, following this recurrence:

where $x_k^u$ means the top bits of slot $k$, $x_{k+1}^l$ means the bottom bits of the next slot, $\mid$ glues them together, and $A$ is a fixed matrix (applied efficiently as a shift plus conditional XOR, shown below).

In plain terms:

1. Take the top bit of the current slot, and the bottom 31 bits of the next slot over. Glue them together into one number, call it $x$.
2. Shift $x$ one bit to the right.
3. If $x$ was an odd number, XOR in a fixed constant, $\texttt{0x9908b0df}$. This is the “twist.”
4. XOR the result with a slot much further ahead, 397 slots over.

Whatever comes out replaces the original slot.

Worked example. Computing the new value for slot 0:

Since $x$ ends in a 0, it’s even, so the twist constant is skipped this time.

That new number, 0x9323659d, isn’t the final output yet, it still needs one more pass.

Temper

The twisted number still has subtle bit patterns in it, so Python runs it through one more fixed sequence of shifts and XORs, called tempering. This is the very last thing the Mersenne Twister does before handing a number back:

Worked example. Running our twisted value through all four passes:

raw output #1 = 0xa3b1799d = 2746317213

That’s a full number out of the Mersenne Twister, start to finish: seed the memory, twist a slot, temper the result. Repeating steps 1c and 1d for the next slot over gives a second number:

raw output #2 = 478163327

Output

Here’s the part that’s specific to Python, not to the Mersenne Twister. A Python float needs 53 bits of precision, but each number out of the twister is only 32 bits. So Python takes two of them and glues them together:

Worked example. Using our two raw numbers from above:

And checking against real Python, seeded identically:

>>> import random
>>> random.seed(42)
>>> random.random()
0.6394267984578837

In conclusion, random.random() gets a seed value, the Mersenne Twister turns it into a stream of 32-bit numbers, and Python glues pairs of them into the float you actually see.

Reference

Below is the twist-and-temper core as C code, matching everything above with the real MT19937 constants:

#define n 624
#define m 397
#define w 32
#define r 31
#define UMASK (0xffffffffUL << r)
#define LMASK (0xffffffffUL >> (w-r))
#define a 0x9908b0dfUL
#define u 11
#define s 7
#define t 15
#define l 18
#define b 0x9d2c5680UL
#define c 0xefc60000UL

uint32_t random_uint32(mt_state* state)
{
    uint32_t* state_array = state->state_array;
    int k = state->state_index;

    int j = k - (n-1);
    if (j < 0) j += n;

    uint32_t x1 = state_array[k] & UMASK;
    uint32_t x2 = state_array[j] & LMASK;
    uint32_t x = x1 | x2;
    uint32_t xA = x >> 1;
    if (x & 1) xA ^= a;

    j = k - (n-m);
    if (j < 0) j += n;

    x = state_array[j] ^ xA;
    state_array[k++] = x;
    if (k >= n) k = 0;
    state->state_index = k;

    uint32_t y = x ^ (x >> u);
             y = y ^ ((y << s) & b);
             y = y ^ ((y << t) & c);
    uint32_t z = y ^ (y >> l);

    return z;
}

Every language’s random module is some version of this same loop, with its own choice of where the seed comes from, and its own formula for turning the raw output into whatever type you asked for.