# Predicting the next pseudorandom value

The archived material for the Stanford University course on cryptography at coursera.org includes a problem where you have to predict the next output of a weak PRG. It can be briefly restated as follows:

\begin{align*} P =&\ 295075153 \\ z_n =&\ x_n \oplus y_n,\qquad \textrm{where}\ 0 \le x_0, y_0 \lt P\ \textrm{(unknown random integers)} \\ x_{n+1} =& (2x_n + 5) \mod P \\ y_{n+1} =& (3y_n + 7) \mod P \\ (z_0, z_1, ..., z_8) =& (210205973, 22795300, 58776750, \\ &\ \ 121262470, 264731963, 140842553, \\ &\ \ 242590528, 195244728, 86752752) \\ z_9 =& ??? \end{align*}

I wrote a small program in C to solve this by brute force, and it finds a solution in just over 2 seconds:

#include <stdio.h>
#define P 295075153L

int main() {
long seq[] = { 210205973, 22795300, 58776750, 121262470, 264731963,
140842553, 242590528, 195244728, 86752752 };
long i, x, y, x0, y0;

for (x0=0; x0<P; x0++) {   /* Try every value of x_0 */
y0 = seq ^ x0;        /* Calculate y_0 from x_0 */
for (x=x0,y=y0,i=1; i<9; i++) {
x = (2 * x + 5) % P;   /* Iterate PRG and check values */
y = (3 * y + 7) % P;
if ((x ^ y) != seq[i]) break;
}
if (i==9) {
printf ("Solution found: x0=%ld, y0=%ld\n",x0, y0);
for (x=x0,y=y0,i=1; i<10; i++) {
x = (2 * x + 5) % P;
y = (3 * y + 7) % P;
}
printf("Next value: %ld\n",x^y);
return 0;
}
}
puts("No solution found");
return 0;
}

I also tried to solve this problem in Python, but ran into two problems:

• Although it does find the same solution, it takes about 3½ minutes to do so. This means it is about a hundred times slower than the C program. Surely this can't be right?

• I'm also not happy with the way I'm checking for a match with the sequence values. In the C code, i is incremented at the end of each for() loop iteration, so when i==9 I can be sure that every value was matched successfully. But this incrementing happens somewhere else in Python, so instead I have to check for i==8 and will report a successful match even if the last comparison fails. How can I do this properly?

from itertools import count

seq = [210205973,22795300,58776750,121262470,264731963,
140842553,242590528,195244728,86752752]
P = 295075153L
for x0 in count():   # better than range(0,P), which is **extremely** slow
if x0==P:
break
y0 = seq ^ x0
x, y = x0, y0
for i in range(1,9):
x = (2 * x + 5) % P
y = (3 * y + 7) % P
if (x ^ y) != seq[i]:
break
if i==8:        # will be true even if seq not matched :-(
break

if i==8:
print "Solution found: x0=%d, y0=%d" % (x0,y0)
x, y = x0, y0
for i in range(1,10):
x = (2 * x + 5) % P
y = (3 * y + 7) % P
print "Next value: %d" % (x ^ y)
else:
print "No solution found"

• A factor of 100 between Python and C is rather common. I recommend you take a look at Cython, capable of getting C speed while still writing Python. – Davidmh Dec 9 '14 at 17:48
• "better than range(0,P), which is extremely slow" In python2 range builds a list, so if P is big it's going to use tons of memory. In python2 you should use xrange (in python3 the old range is done and xrange was renamed range). – Bakuriu Dec 9 '14 at 21:05
• BTW I believe the problem can be solved with just paper and pencil. There are a few theorems out there about equalities modulo a prime P, so I'm pretty sure we can do much better than brute-forcing. Once you find a better algorithm the difference in time probably will be much smaller for the simple fact that the solution will be found in a matter of instants in both programs. – Bakuriu Dec 9 '14 at 21:09

## 3 Answers

Here's a refactor of mjolka's version:

def generate(x, y, P, seq):
for i in seq:
x = (2 * x + 5) % P
y = (3 * y + 7) % P
if x ^ y != i:
return -1

x = (2 * x + 5) % P
y = (3 * y + 7) % P

return x ^ y

def solve():
seq0 = 210205973
seq = [22795300,  58776750,  121262470, 264731963,
140842553, 242590528, 195244728, 86752752]

P = 295075153
for x0 in range(0, P):
y0 = seq0 ^ x0

if generate(x0, y0, P, seq) != -1:
print("Solution found: x0={}, y0={}".format(x0, y0))
value = generate(x0, y0, P, seq)
print("Next value: {}".format(value))
return

print("No solution found")

solve()


This solves your second problem by having a function "extract" the control flow, and use a sentinel as the return to signal failure to match.

I used -1 as the "not found" return value because it's significantly faster than the cleaner alternative of returning None on PyPy.

Times:

$time pypy p.py Solution found: x0=89059908, y0=164204369 Next value: 86752752 pypy p.py 1.86s user 0.01s system 99% cpu 1.874 total$ time pypy3 p.py
Solution found: x0=89059908, y0=164204369
Next value: 86752752
pypy3 p.py  2.00s user 0.00s system 99% cpu 2.012 total


This is actually slightly faster than mjolka's:

$time pypy p.py Solution found: x0=89059908, y0=164204369 Next value: 231886864 pypy p.py 2.12s user 0.01s system 99% cpu 2.136 total  This was significantly slower than the C compiled with -O3: $ gcc -O3 c.c -o c

$time ./c Solution found: x0=89059908, y0=164204369 Next value: 231886864 ./c 0.28s user 0.00s system 99% cpu 0.282 total  but if the value of P is instead decided at runtime, the speeds are very close: $ time echo 295075153 | ./c
Solution found: x0=89059908, y0=164204369
Next value: 231886864
echo 295075153  0.00s user 0.00s system 0% cpu 0.001 total
./c  1.67s user 0.00s system 99% cpu 1.671 total


(1.67s vs 1.87s). I used C++ to get input, because I'm lazy:

  std::cin >> P;// = 295075153L;


If you accept an uglier form of looping, this is as fast, removes the duplication of x = ...; y = ... and doesn't split seq:

def generate(x, y, P, seq):
for i in range(1, 10):
x = (2 * x + 5) % P
y = (3 * y + 7) % P

if i == 9:
return x ^ y

if x ^ y != seq[i]:
return -1

def solve():
seq = [210205973, 22795300,  58776750,  121262470, 264731963,
140842553, 242590528, 195244728, 86752752]

P = 295075153
for x0 in range(0, P):
y0 = seq ^ x0

if generate(x0, y0, P, seq) != -1:
print("Solution found: x0={}, y0={}".format(x0, y0))
value = generate(x0, y0, P, seq)
print("Next value: {}".format(value))
return

print("No solution found")

solve()


Times:

$time pypy p.py Solution found: x0=89059908, y0=164204369 Next value: 231886864 pypy p.py 1.89s user 0.01s system 99% cpu 1.905 total$ time pypy3 p.py
Solution found: x0=89059908, y0=164204369
Next value: 231886864
pypy3 p.py  2.00s user 0.00s system 99% cpu 2.003 total

• I don't like that you have to split seq and that you call x = (2 * x + 5) % P and the other twice. – njzk2 Dec 9 '14 at 18:59
• @njzk2 I split seq because passing seq[1:] lots of times proved slow. I personally also find applying x = ...; y = ... simpler than the other techniques despite duplication, but I've added another option that avoids that anyway (see update). If you're willing to bear the small runtime cost, you could instead just make a function to do an iteration and call that. – Veedrac Dec 9 '14 at 19:32
• you can test z = x^y before updating x and y, removing the need to slice seq. that makes one useless test, though. – njzk2 Dec 9 '14 at 19:34
• @njzk2 Since almost all attempts succeed 0 times, adding one extra test like that actually costs far more than you'd expect. I think it also hurts unrolling. The code took about 3 times as long when I tried that. – Veedrac Dec 9 '14 at 19:53
• good point. I may try with an iterator, then (to consume the first item as seq0 and loop on the rest) – njzk2 Dec 9 '14 at 20:01

To address your question:

In the C code, i is incremented at the end of each for() loop iteration, so when i==9 I can be sure that every value was matched successfully. But this incrementing happens somewhere else in Python, so instead I have to check for i==8 and will report a successful match even if the last comparison fails. How can I do this properly?

In Python a for loop can have an else clause that is executed when the loop ends normally -- not via break. In your case you could simply change if i==8: into else:. The break under else: will then break the outer loop.

There are a few code issues that tools like flake will pick up on, but I'll just focus on performance.

First up, switching to pypy reduces the runtime on my machine from ~3 minutes to 16 seconds.

⚡ time pypy orig.py
Solution found: x0=89059908, y0=164204369
Next value: 231886864
pypy orig.py  15.79s user 0.07s system 99% cpu 15.915 total


Removing the L gives us another drastic speedup:

⚡ diff orig.py orig2.py
5c5
< P = 295075153L
---
> P = 295075153

⚡ time pypy orig2.py
Solution found: x0=89059908, y0=164204369
Next value: 231886864
pypy orig2.py  3.93s user 0.04s system 99% cpu 3.988 total


Changing the control flow gives us a smaller boost

def solve():
seq = [210205973, 22795300, 58776750,
121262470, 264731963, 140842553,
242590528, 195244728, 86752752]
P = 295075153
for x0 in xrange(0, P):
y0 = seq ^ x0
x, y = x0, y0
for i in xrange(1, 9):
x = (2 * x + 5) % P
y = (3 * y + 7) % P
if (x ^ y) != seq[i]:
break
if i == 8:
print "Solution found: x0=%d, y0=%d" % (x0, y0)
x, y = x0, y0
for i in xrange(1, 10):
x = (2 * x + 5) % P
y = (3 * y + 7) % P
print "Next value: %d" % (x ^ y)
return
print "No solution found"

solve()

⚡ time pypy orig3.py
Solution found: x0=89059908, y0=164204369
Next value: 231886864
pypy solve.py  2.60s user 0.03s system 99% cpu 2.633 total


Which gets us very close to the C solution

⚡ time ./a.out
Solution found: x0=89059908, y0=164204369
Next value: 231886864
./a.out  2.36s user 0.01s system 99% cpu 2.369 total


Forgot to compile with optimisations, which puts the Python version ~6x slower than C:

⚡ time ./a.out
Solution found: x0=89059908, y0=164204369
Next value: 231886864
./a.out  0.44s user 0.00s system 99% cpu 0.448 total