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 nth output is an z reduced mod m. So our task is to solve
x = an 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 = loga (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 = 231 – 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)
2 thoughts on “Predicting when an RNG will output a given value”
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.