\(\mathbb{F}_2\)-Linear PRNGs
See Also
An operation \(f(x)\) can be described as "\(\mathbb{F}_2\)-linear" if, when treating the bits of \(x\) as a binary vector (one whose elements are of the finite field \(\mathbb{F}_2\)), \(f(x)\) can be written as \(f(x) = xA\) for some matrix \(A\) (alternatively, \(f(x)^T = A^Tx^T\)). A PRNG is described as \(\mathbb{F}_2\)-linear if its state transition function has this property.
Examples of \(\mathbb{F}_2\)-Linear Operations
Over \(\mathbb{F}_2\), addition becomes equivalent to binary XOR and multiplication becomes equivalent to binary AND. This results in many bitwise operations being \(\mathbb{F}_2\)-linear. Due to the properties of matrix multiplication, each row \(i\) of \(A\) can be thought of as encoding what the \(i\)th bit of \(x\) contributes to \(f(x)\). The properties of matrix multiplication also imply composed \(\mathbb{F}_2\)-linear operations are also \(\mathbb{F}_2\)-linear (\(xA_1 A_2 A_3... = xB\)).
For clarification, binary vectors are constructed from the least significant to the most significant bit, that is, \(x_0 =\) x & 1, \(x_1 =\) (x >> 1) & 1, etc..
Identity
The identity operation \(f(x) = x\) can be represented as the identity matrix \(I\). This is trivial to show given that \(xI = x\). This is also shown by the fact that each row \(i\) of \(I\) contains only one non-zero entry at column \(i\). The contribution to each bit \(i\) of \(f(x)\), then, is only made up of that same bit \(i\) of \(x\).
Masking (AND by constant)
Masking by a constant \(f(x) =\) x & c is the same as setting a constant set of bits \(S\) to \(0\) and leaving the rest as is. \(A\), then, should look like the identity matrix for rows \(i \notin S\) and \(0\) for rows \(i \in S\). This can be constructed by placing the bits of c along the diagonal of a square zero matrix.
For example, masking a 4-bit \(x\) by 0xb = 0b1011 is represented by the matrix
XOR by a variable
XORing two functions of \(x\) \(f(x) = g(x)\) ^ \(h(x)\) is the same as adding their \(\mathbb{F}_2\) representations \(Fx + Gx = (F + G)x\). The matrix representing \(f(x)\), then, is simply \(A = F + G\). For \(f(x) = x\) ^ \(g(x)\) this is simply \(A = I + G\).
Bitwise shift
Shifting by a constant \(f(x) =\) x >> a (where \(a < 0\) implies a left shift by \(-a\)) is the same as setting each bit of the output \(f(x)_i\) to a corresponding bit of the input \(a\) away \(x_{i+a}\) (if it exists, otherwise it is set to \(0\)). Each row \(i\) of \(A\) should reflect that the only contribution is from bit \(i+a\) meaning each row consists of only one non-zero entry at column \(i+a\).
Notably, given a matrix \(R\) representing x >> 1, x >> a can be constructed by composing the operation \(A = R^a\). The same goes for a left shift matrix \(L\) where x << a is represented by \(L^a\).
For example, right shifting a 4-bit \(x\) by 1 is represented by the matrix
and left shifting by
It is easy to see visually that these can also be obtained by shifting the columns of the identity matrix.
Bitwise rotate
Bitwise rotation by a constant \(f(x) =\) rotl(x,a) (where \(a < 0\) implies a right rotation by \(-a\)) can be represented by first decomposing the rotation into shifts rotl(x, a) = (x << a) | (x >> (d - a)) where \(d\) is the bitlength of \(x\). By the nature of rotl, the shifts result in no overlapping bits, so the OR operation is equivalent to either addition or XOR. This gives rotl(x, a) = lshift(x, a) ^ rshift(x, d-a), or, \(A = L^a + R^{d-a}\). A similar construction can be found by noticing the operation maps \(x_{i}\) to \(f(x)_{(i+a) \bmod d}\).
Notably, given a matrix \(R_R\) representing rotl(x, -1), rotl(x, -a) can be constructed by composing the operation \(A = R_R^a\). The same goes for a left rotation matrix \(L_R\) where rotl(x, a) is represented by \(L_R^a\).
For example, rotating a 4-bit \(x\) left by 1 is given by
and right rotating by
It is easy to see visually that these can also be obtained by rotating the columns of the identity matrix.
Xorshift
Xorshift by a constant \(f(x) =\) x^(x >> a) (where \(a < 0\) implies a left shift by \(-a\)) is simply the XOR of the identity by a left or right shift. This can be represented by simply adding the identity matrix to a shift matrix \(X_a = I + R^a\) or \(X_a = I + L^{-a}\).
Conditional XOR by a constant
Conditial XOR by a constant, that is, XORing by some constant \(c\) if some bit \(k\) of \(x\) is set \(f(x) = x\) ^ (((x >> k) & 1) * c) or
if ((x >> k) & 1) {
x ^= c;
}
can be represented by first representing the conditional (x >> k) & 1 ? c : 0. \(k\) is the only bit that contributes to the output, and it contributes precisely the bits of \(c\) when its set. The matrix \(C\), then, should have \(0\)s everywhere but row \(k\) which should contain the bits of \(c\). The full operation can be constructed by XORing the identity operation by this conditional \(f(x) = x\) ^ \(g(x)\) which gives \(A = I + C\).
For example, XORing a 4-bit \(x\) by 0xb = 0b1011 when \(x_0\) is set \(f(x) = x\) ^ ((x & 1) * 0xb) is represented as follows
Direct Construction
If an operation is known to be \(\mathbb{F}_2\)-linear, its corresponding matrix can be constructed by directly testing the contribution from each input bit. If the input \(x\) consists of only \(1\) set bit at index \(i\), multiplication by \(A\) will result in an output of \(A\)'s \(i\)th row, or, \(A_i = f(2^i)\).
Efficient Jumping
Efficiently jumping an \(\mathbb{F}_2\)-linear PRNG requires finding an easy to compute function \(f(x, n) = xA^{n}\) that is equivalent to \(n\) applications of the state transition function. \(A^n\) can be calculated efficiently via exponentiation by squares. This achieves time complexity \(O(\log n)\) but requires expensive matrix computations.
A more efficient method can by achieved by utilizing the Cayley-Hamilton Theorem which states that a square matrix \(A\) satisfies its own characteristic equation \(p(A) = 0\) where \(p\) is the characteristic polynomial of the matrix. Define a "jump polynomial" \(j_n(x) = x^n\) which can be used to write the goal function \(f(x, n) = xj_n(A)\). \(j_n(x)\) can be decomposed into the form \(q(x)p(x) + r(x)\) for some polynomials \(q(x)\) and \(r(x)\). Evaluating this at \(A\) results in \(q(A)p(A) + r(A) = r(A)\). This reveals that any jump polynomial can be reduced by the characteristic polynomial and still produce the same matrix \(x(j_n(z) \bmod p(z))(A) = xA^{n}\). A reduced jump polynomial will have degree \(\deg(p) - 1\) and can be evaluated with that many state transitions via Horner's Method. Each step requires only multiplication by \(A\) (which can be fully replaced with the actual state transition function to avoid matrices) and addition between state vectors (which becomes XOR).
For a constant \(n\), \(j_n\) can be precomputed and the time complexity of jumping is \(O(\deg(p))\) which is constant for a chosen PRNG. If \(n\) is not constant, computing it requires \(\mathbb{F}_2\) polynomial exponentiation which is much faster than matrices but still \(O(\log n)\).
Constant time can be achieved for variable \(n\) by utilizing the fact that jumps can be decomposed such that \(f(x, n+m) = f(f(x, n), m)\). For a PRNG with period \(P\), any \(n\) can be reduced \(\bmod P\) so only \(n \in [0, P - 1]\) need to be considered. If \(n\) can be decomposed into a sum \(\sum_{i=1}^{N}n_i\) where each \(j_{n_i}\) is known, the jump can be performed in \(N(\deg(p) - 1)\) state transitions via repeated applications of the corresponding jump polynomials.
An easy to calculate and work with decomposition of \(n\) is given by its base-2 representation \(\sum_{i=1}^{N}2^{i-1}c_i\) where \(N = \lceil \log_2 P \rceil\). This results in \(\lceil \log_2 P \rceil(\deg(p) - 1)\) state transitions which is constant time and does not require matrix or polynomial computations.
An example for a 64-bit PRNG could be implemented as follows:
uint64_t JUMP_POLYNOMIALS[64] = {...}; // J[i] jumps 2^i times
void jump(uint64_t &state, uint64_t n) {
for (int i = 0; n; n >>= 1, i++) {
if (n & 1) {
// jump 2^i
uint64_t result = 0;
uint64_t J = JUMP_POLYNOMIALS[i];
// evaluate via Horner's method
for (; J; J >>= 1) {
if (J & 1) {
result ^= state;
}
next_state(state);
}
state = result;
}
}
}
State Recovery
If the PRNGs output function (or a part of it) is also \(\mathbb{F}_2\)-linear, the linearity can be leveraged to recover the internal state from measured outputs.
If the \(i\)th output can be written as \(xA^iO_i\) for some output matrix \(O_i\), the state transition matrix \(A\), and the initial state \(x\), start by incorporating the state transition matrix into the output \(O_i := A^iO_i\). These matrices can be concatenated horizontally to produce one large matrix \(O\) that maps an initial state to a vector containing the bits of all the outputs. Once sufficient outputs are collected such that the system is fully determined, a right inverse \(O^{-1}\) exists such that \(xOO^-1 = x\). As \(xO\) is exactly the observed outputs and \(O\) is easily constructed, the internal state \(x\) can be computed by simply inverting the matrix.