Arithmetic background of this RNG code

Question

I am trying to figure out the mathematical background of the random number generation of an old video game. It does iterations where it updates a 33-bit state consisting of the variables z (32-bit) and some flag f (1-bit). It generates a 'random' 8-bit value per iteration by just taking the lowest byte of z. This is the algorithm in C++ style pseudocode:

    unsigned int z = SOME_SEED; // state variable
    unsigned int f = 0;         // LSB of previous state z
    unsigned int t;             // temp variable

    while(someRunningCondition)
    {
        t = ((z >> 1) + (f << 31)) ^ (z << 12);
        f = z & 1;
        z = t ^ (t >> 20);

        outputRandomValue(z & 0xFF);
    }

The iteration visits every state other than z = 0, f = 0 and thus has a period of $2^{33}-1$ and the random samples have uniform distribution. The code looks similar to efficient implementations of finite field multiplication in cryptography, but I am too inexperienced in the topic to figure it out myself. If someone could break down the math underlying this code for me, that would be great.

Answer 1

Whether 2ⁿ − 1 is prime isn’t very important if we look for multiplication by certain constant in F_2ⁿ . Multiplicative groups of finite fields are always cyclic, and any cyclic group has elements of the full order, prime or composite. As follows from inspection the code, mathematics of this example is iteration of a linear operator on 33-dimensional vector space over F₂ . The only dubious fragment (z >> 1) + (f << 31) reveals addition modulo 2 as well, because the former term has most significant bit zero, whereas the latter has zeros anywhere but in the MSB. As any field extension, F_2³³ is isomorphic to a vector space, but to find a multiplication in this operator we must specify an isomorphism, i.e. specify how is F_2³³ encoded with these z and f.

We demonstrated it encode a multiplication in F_2³³ as follows. First, one observation:

t = … (z >> 1) … ^ (z << 12); relative offset is 13 positions.

z = t ^ (t >> 20); relative offset is 20 positions.

But 13 + 20 = 33. There is well-known representation of F_pⁿ (p = 2, n = 33 in our case) as a quotient ring of polynomials over F_p over the ideal generated by some irreducible polynomial of degree n; we need to guess the ideal and the factor (multiplication to which does the operator represent).

Matrix of the operator that transforms f and z on each iteration in coordinates (f, z_MSB, … , z_LSB) is:

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1
1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0

Bands of ones progressing right-downwards are evident. We also see, when such band wraps over the right boundary, it starts

1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 …

on the left side, and nothing more is created there. Transpose of this matrix is multiplication by X³² in the quotient ring F₂[X]/〈X³³ − X¹³ − 1〉 in the standard basis (from 1 to X³²).

Added on original poster’s insistence: Why is F₂[X]/〈X³³ − X¹³ − 1〉 isomorphic to F_2³³? Because we can prove numerically that order of multiplication-by-X³² operator is 2³³ − 1, hence any non-zero element of the quotient ring is a power of X³². It implies that the ring is a field and X³² generates its multiplicative group. There is only one field, up to isomorphism, with given finite number of elements.

Namely:

 X⁰ × X³² = X³² (1st row shown above; corresponds to f = z & 1)
 X¹ × X³² ≡ X¹³ + X⁰ (2nd row shown above; expression for MSB of z)
 X² × X³² ≡ X¹⁴ + X¹ (3rd row shown above)
   ⋮
 X²⁰ × X³² ≡ X³² + X¹⁹
 X²¹ × X³² ≡ X²⁰ + X¹³ + X⁰
   ⋮

Then, what means the transpose clause? It means that elements of F_2³³ are not 33-bit (f, z) vectors (the state), but F₂-linear functionals on them. For example, z & 0xFF can be thought of as eight bit functionals on the state. Multiplication by X³² iterates these functionals, that are then applied to initial seed. Behaves similarly stochastic; a mathematical difference there is, but a practical difference there is not.

Does it mean this operator isn’t a direct representation of some multiplication in F_2³³ ? Of course, not. To obtain a matrix with the same structure of columns (instead of rows) we can just permute the basis. Namely, let’s write it in coordinates starting from 12th bit of z (the coefficient at 2²⁰ in unsigned int), then counting 11th, 10th, 9th… bits of z up to MSB, then f, then 32th bit of z, and so on. Then, we see exactly transpose of the original matrix. IMHO all such transposes of multiplications, in all finite fields, are “the same” (i.e. same conjugacy class in GL_n(F_p)) as operators of (possibly other) multiplications in the same field, but I’m reluctant to investigate.