Breaking Linear Congruential Generators

There is a nice trick to solve for $m.$ Let $Y_n = X_n-X_{n-1}.$ Then $Y_n \equiv (aX_{n-1}+b) - (aX_{n-2}+b) \equiv a(X_{n-1}-X_{n-2}) \equiv aY_{n-1} \pmod{m}.$ Now let $Z_n = Y_nY_{n+2}-Y_{n+1}^2.$ Then $Z_n$ is congruent to $0$ mod $m.$

(Proof: each $Y_k$ is congruent to $a Y_{k-1},$ so we can express everything in terms of $Y_n:$ $Z_n = Y_nY_{n+2} - Y_{n+1}^2 \equiv Y_n (a^2 Y_n) - (aY_n)(aY_n) \equiv 0 \pmod m.)$

Since each $Z_n$ is a multiple of $m,$ the idea is to generate enough values of $Z_n$ and take their gcd, and expect that the answer equals $m.$ (With sufficiently many $Z_n,$ we expect there to be no larger common factor than the one that is guaranteed, namely $m.$ )

Anyway, the gcd of $Z_1,\ldots,Z_6$ is $m=1000000007.$

Given $m,$ it is easy enough to solve for $a$ and $b$ ; just use $\begin{aligned} aX_1+b &\equiv X_2 \\ aX_2+b &\equiv X_3 \end{aligned}$ to get two equations in the unknowns $a,b,$ whose solutions are $a \equiv \frac{X_2-X_3}{X_1-X_2}, \quad b \equiv X_2-aX_1.$ (These computations are all done mod $m.$ When I did them, I cheated and used that $m$ was prime, which let me compute $1/d$ as $d^{m-2} \pmod m.$ Without the cheat, I would have had to write an extended Euclidean algorithm function to compute the modular inverse of $X_1-X_2$ .)

This gives $a=1664525$ and $b=13904216,$ whence the next term is $1664525 \cdot 200183345 + 13904216 \equiv \fbox{193907871} \pmod{1000000007}.$

This solution written in python 3 based on the same idea of Patrick Corn solution.

Extended Euclidean algorithm is from wikibooks

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

from fractions import gcd

Xn = [720555190, 133143292, 350469176, 715002068, 822810950, 400865843, 226553034, 200183345]

def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception('modular inverse does not exist')
    else:
        return x % m

def solve(Xn):
    "Breaking Linear Congruential Generators"
    Yn = [Xn[k+1]-Xn[k] for k in range(len(Xn)-1)]

    Zn = [abs(Yn[k+2]*Yn[k] - Yn[k+1]*Yn[k+1]) for k in range(len(Yn)-2)]

    m = Zn[0]
    for x in Zn[1:]:
        m =  gcd(m, x)

    a = (modinv(Yn[0] + m, m)*Yn[1])%m
    b = (Xn[1] - a *  Xn[0])%m

    x = Xn[0]
    for i in range(8):
        if x != Xn[i]:
            raise Exception('Failed to solve')
        x = (a*x + b)%m
        return x

print(solve(Xn))