A few days ago I wrote about how to pick the seed of a simple random number generator so that a desired output came *n* values later. The number *n* was fixed and we varied the seed. In this post, the seed will be fixed and we’ll solve for *n*. In other words, we ask when a pseudorandom sequence will produce a given value.

In principle you could just run the RNG until you get the output you’re looking for, but we’ll assume such a brute force approach is not feasible or at least not fast enough.

If a LCG (linear congruential generator) has seed *z*, multiplier *a*, and modulus *m*, then the *n*th output is a^{n} *z* reduced mod *m*. So our task is to solve

*x* = *a*^{n} *z* mod *m*

for *n*. If we forget for a moment that we’re working with integers mod *m*, we see that the solution is

*n* = log_{a} (*x* / *z*)

We can actually make this work if we interpret division by *z* to mean multiplication by the inverse of *z* mod *m* and if we interpret the logarithm to be a discrete logarithm. For more on discrete logarithms and one algorithm for computing them, see this post.

In an earlier post I used *a* = 742938285 and *m* = 2^{31} – 1 = 2147483647. We set *n* = 100 and solved for *z* to make the 100th output equal to 20170816, the date of that post. It turned out that *z* = 1898888478.

Now let’s set the seed *z* = 1898888478 and ask when the LCG sequence will take on the value *x* = 20170816. Of course we know that *n* will turn out to be 100, but let’s pretend we don’t know that. We’ll write a little Python script to find *n*.

I expect there’s a simple way to compute modular inverses using SymPy, but I haven’t found it, so I used some code from StackOverflow.

The following Python code produces *n* = 100, as expected.

from sympy.ntheory import discrete_log 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 a = 742938285 z = 1898888478 m = 2**31 - 1 x = 20170816 zinv = modinv(z, m) n = discrete_log(m, x*zinv, a) print(n)

Sympy has ‘gcdex’, ‘half_gcdex’, ‘invert’ and ‘sympy.core.numbers.mod_inverse’ and you can also use

zinv = (FF(m)(1) / FF(m)(z)).to_int()

SymPy has the mod_inverse() function, which is the same as your modinv().

Of course, it would be cool if you could just use solveset(x – Mod(a**n*z, m), n, S.Integers), but it isn’t implemented yet.