Attack of the Middle Square Weyl Sequence PRNG

Question

The Middle Square Weyl Sequence PRNG was proposed as illustration of properties of CSPRNGs in this answer (now deleted). This PRNG was introduced by Bernard Widynski in Middle Square Weyl Sequence RNG (arXiv, 2017). That was without claim of cryptographic security, but v5 of the paper conjectures that it would "be difficult to determine the internal state by examining the outputs".

The PRNG has a 128-bit state consisting of two 64-bit variables $x$ and $w$. The PRNG input consists of the initial $x_0$ and $w_0$, and an odd value $s$. The state's evolution and 32-bit pseudorandom outputs $r_i$ are defined (for $i\ge0$) by $$\begin{align} r_i&\gets x_i\bmod2^{32}\\ w_{i+1}&\gets(w_i+s)\bmod 2^{64}\\ x_{i+1}&\gets((x_i^2+w_{i+1})\bmod 2^{64})\ggg32 \end{align}$$ where $\ggg32$ swaps the halves of its 64-bit left argument.

^{[The notation was modified so that $r_0$ is output, for consistency with an answer]}

The paper gives a reference C99 implementation (with input $x_0=w_0=0$ and a certain $s$; the output starts at $r_1$)

uint64_t x = 0, w = 0, s = 0xb5ad4eceda1ce2a9;
inline static uint32_t msws() {
   x *= x; x += (w += s); return x = (x>>32) | (x<<32);
}

Assume we have a $r_i$ for $0\le i\le9$. How can the rest of the sequence be predicted efficiently? That would be a total break.

If not feasible, lesser breaks (only distinguishing the output from random, or requiring more plaintext, or working only for a fraction of the inputs) are also welcome.

Answer 1

Lets look at the computations that the function does to the internal state.

We'll denote the state of the variables when we start the attack as $x_0 = (a, b)$ and $w_0 = (c, d)$, where I will use the notation $(i, j)$ to mean the 64 bit number which has $i$ as the upper 32 bits, and $j$ as the lower 32 bits (and so $a, b, c, d$ are 32 bits each). Note that $x_0, w_0$ need not be the initial state of the PRNG, instead it's the state when the initial output was generated.

Then, the previous step exposes $b$ (and so that is known to the attacker).

Then, it squares $x_0$ modulo $2^{64}$, replacing $x_0 = (a, b)$ with $x_0^2 = (2 a b + (b^2 >> 32), b^2\ \bmod 2^{32})$. As the attacker knows $b$, the attacker can replace this with the linear operation $x_0^2 = (h_0 a + h_2, h_1)$, where $h_0, h_1, h_2$ are known (to the attacker) constants (as it gives the same mapping from previous state to new state).

Then, the algorithm adds $w$ to $x$, that can be modeled as:

$x_0^2 + w_0 = (h_0 a + h_2 + c + \epsilon_0, h_1 + d)$, where $\epsilon_0$ is either 0 or 1, depending on whether a carry from the lower 32 bits happened.

Then, it adds the constant $s$ to $w$, this can be modeled as $w + s = (c + s_{hi} + \epsilon_1, d + s_{lo})$, where $\epsilon_1$ is either 0 or 1, depending on whether a carry from the lower 32 bits happened (and $s = (s_{hi}, s_{lo})$

Then, it rotates $x$, giving $x_1 = (h_1 + d, h_0a + h_2 + c + \epsilon_0)$, and outputs the lower 32 bits $h_0a + h_2 + c + \epsilon_0$.

The new state is $x_1 = (h_1 + d, h_0a + h_2 + c + \epsilon_0)$ and $w_1 = (c + s_{hi} + \epsilon_1, d + s_{lo})$; that is, it is one of 4 linear functions of the previous state (depending on the values of $\epsilon_0, \epsilon_1$)

Continuing on for three more iterations, the resulting states will remain one of a small number of known linear functions of the initial state; the output will also be a linear function of the state. With four outputs, the attacker can iterate through the possible $\epsilon$ values, and solve for the initial states. If you go through the equations, you'll find that you can recover $c, d$ from the third and fourth outputs; using the value $c$ gives you the value of $a$ from the second output (and you're given $b$ as the first output).