prg_cracker.py

#!/usr/bin/python2.7
# -*- encoding: utf-8 -*-
from __future__ import division

import sys

P = 0xFFFFFFFFFFFFFFC5 # The largest 64-bit prime

def f(n):
    return (2*n) ^ ((3*n + 7) % P)

# A function related to f that is easy to invert
def _f(n):
    return (2*n) ^ (3*n + 7)

# Inverse of _f
def _g(n):
    i, r = 1, 0
    while i <= n:
        i, r = i * 2, r + ( ((n ^ _f(r + i)) & i) ^ i )
    return r

def f_inverse(n):
    """Try reasonably hard to find a value y such that f(y) == n.
    
    The key observation is that f(y) and (_f(y) % P) usually
    differ only in their low-order bits, hence _g( f(y) + kP ) is
    usually close to y for some small k.
    """
    for k in (0,1,2,-1,3,4,5,6,8): # Values determined experimentally
        if n + k*P < 0: continue
        x = _g(n + k*P) % P
        for i in xrange(0x10000):
            y = x ^ i
            if f(y) == n:
                return y
    return None

def run_prng(x, y, init_i, n):
    for i in range(init_i, init_i + n):
        print "Output #{}: {}".format(i+1, x^y)
        x, y = (2*x + 5) % P, (3*y + 7) % P

# These values were also determined experimentally. The ones of the form
# P + x will depend on the value of P, I assume. The code to find these
# values is in ts.py in this gist.
ts = (5, 13, 29, 61, 125, P + 123, P + 251, P + 507, P + 1019, P + 2043)

def crack(ns):
    for i, (a, b) in enumerate(zip(ns, ns[1:])):
        for t in ts:
            f_y = (2*a) ^ b ^ t
            y = f_inverse(f_y)
            if y is None:
                continue
            
            x = a^y
            if ((2*x + 5) % P) ^ ((3*y + 7) % P) == b:
                run_prng(x, y, i, 11)
                return

ns = map(int, sys.argv[1:])
crack(ns)

ts.py

import os, struct

P = 0xFFFFFFFFFFFFFFC5 # The largest 64-bit prime

def rand():
    while True:
        x, = struct.unpack('Q', os.urandom(8))
        if x < P:
            return x

p = {}
for i in range(100000):
    x = rand()
    n = (2*x) ^ ((2*x + 5) % P)
    p[n] = p.get(n, 0) + 1

s = set()
for n, occurences in p.iteritems():
    if occurences > 1000:
        s.add(n)

print "(" + ", ".join([ str(n) if n<P else "P + " + str(n-P) for n in sorted(s) ]) + ")"