LGC

Linear Congruential Generator.

n (modulus), a (multiplier), b (increment), Seed X0

Xi+1 = aXi + b (mod n)

Attack with n, a known

b = Xi+1 - aXi (mod n)

Attack with n known

Xi+1 = aXi + b (mod n)
Xi+2 = aXi+1 + b (mod n)

a = (Xi+2 - Xi+1) * (Xi+1 - Xi)^-1 (mod n)

Demonstration

Xi+1 = aXi + b (mod n)
Xi+2 = aXi+1 + b (mod n)

Xi+2 - Xi+1 = aXi+1 +b - (aXi + b) (mod n)
Xi+2 - Xi+1 = aXi+1 +b -aXi -b (mod n)
Xi+2 - Xi+1 = a(Xi+1 - Xi) (mod n)
(Xi+2 - Xi+1) * (Xi+1 - Xi)^-1 = a(Xi+1 - Xi)*(Xi+1 - Xi)^-1 (mod n)

a = (Xi+2 - Xi+1) * (Xi+1 - Xi)^-1 (mod n)

Attack with nothing known

Take the differences 
Tn = Xn+1 - Xn                 For every n we know - 1

Take more sequence
Un = | Tn+2 * Tn - Tn+1^2 |    For every n we know - 3

With high probability we have that
n = gcd(U1, U2, U3, ...)       For each U we have

Demonstration

Tn = Xn+1 - Xn
   = (aXn + b) - (aXn-1 + b) mod n
   = aXn + b - aXn-1 - b mod n
   = aXn - aXn-1 mod n
   = a(Xn - Xn-1) mod n
   = aTn-1 mod n

Tn+2 = aTn-1 mod n
     = a^2Tn mod n

Un = | Tn+2 * Tn - Tn+1^2 | 
   = | a^2Tn * Tn - (aTn)^2 | mod n
   = | a^2Tn^2 - (aTn)^2 | mod n
   = | (aTn)^2 - (aTn)^2 | mod n
   = 0 mod n
   That is, multiple of n

Getting different multiples of n

Now, the gcd of two random multiples of n will be n with probability approximately 0.61. The more U numbers you have, the higher the probability (usually 4/5 U are enough).

Code

import math
from functools import reduce

def compute_next(output, modulus, a, b):
    return  ((output*a) + b) % modulus
def compute_increment(outputs, modulus, a):
    b = (outputs[1] - outputs[0] * a) % modulus
    return b
def compute_multiplier(outputs, modulus):
    term_1 = outputs[1] - outputs[2]
    term_2 = pow(outputs[0] - outputs[1], -1, modulus) 
    a = (term_1 * term_2) % modulus
    return a
def compute_modulus(outputs):
    ts = []
    for i in range(0, len(outputs) - 1):
        ts.append(outputs[i+1] - outputs[i])
        
    us = []
    for i in range(0, len(ts)-2):
        us.append(abs(ts[i+2]*ts[i] - ts[i+1]**2))

    modulus =  reduce(math.gcd, us) 
    return modulus
    
n = compute_modulus(outputs)
a = compute_multiplier(outputs, n)
b = compute_increment(outputs, n, a)