Predicting the next pseudorandom value

Question

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[0] ^ 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[0] ^ 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[8] 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"

Answer 1

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[0] ^ 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

Answer 2

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.

Answer 3

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[0] ^ 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