# Any speed up methods for finding the index of a random value produced by the Inversive congruential generator?

The Inversive congruential generator produces random values with: $$x_{n+1} = a\cdot x_{n}^{-1} + b \mod P$$

(special case if $$x_n=0$$ -> $$x_{n+1}=b$$)
starting with an initial value $$x_0$$

With well choosen $$a,b,P$$ it produces all values in $$\mathbb{F}_P$$ (which are $$[0..P-1]$$)
That's the case if there is no r with $$r^2 \equiv 4 a + b^2 \mod P$$ with $$P$$ a prime.

Question: Is there any way to compute the index $$n$$ for a given value $$v \in \mathbb{F}_P$$ with $$x_n=v$$ ? (with $$a,b,P$$ well chosen and known)
Or is there any way to compute it faster than brute force?

For the $$n$$'th random value a equation exists with: $$x_n=(q_{n+1}\cdot x_0 + a\cdot q_n)(q_{n}\cdot x_0 + a\cdot q_{n-1})^{-1} \mod P$$ with $$q$$ a sequence again: $$q_0=0$$ $$q_1=1$$ $$q_{n}=a\cdot q_{n-2}+b\cdot q_{n-1} (\mod P)$$

q is similar to the Fibonacci sequence. Without modulo there is an equation for the $$n$$'th element as well: $$q_n= \frac{ (\frac{b +\sqrt{4 a + b^2} }{2})^n - (\frac{b - \sqrt{4 a + b^2}}{2})^n}{\sqrt{4 a + b^2}}$$ (I don't know if there also exists a form modulo $$P$$)

Update: $$q_n$$ equation modulo P
with:
$$r^2 \equiv 4 \cdot a+b^2 \mod P$$
$$1 \equiv t \cdot 2 \mod P$$

$$q_n = ((b+r)^n-(b-r)^n)\cdot t^n \cdot r^{-1} \mod P$$ This equation also works in $$\mathbb{F}_P$$

However for the target values $$a,b$$ there exists no such root $$r$$.
(For other $$a,b$$ which can not produce all values it does work)

Update 2: simplification
For simplification we can assume $$x_0=0$$. With this the equation for $$x_n$$ would be: $$x_n=(a\cdot q_n)(a\cdot q_{n-1})^{-1} \mod P$$ $$x_n=q_n\cdot q_{n-1}^{-1} \mod P$$

Any way to compute index $$n$$ for given value $$v=x_n$$?

Update 3: some test in wolframalpa
wolfram alpha has a solution for the non-modulo version (with $$x0=0$$): https://www.wolframalpha.com/input/?i=solve+(((b%2Br)%2F2)^n-((b-r)%2F2)^n)%2F(((b%2Br)%2F2)^(n-1)-((b-r)%2F2)^(n-1))%3Dv+for+n?

$$n = \frac{ \log(\frac{(b - r) (b + r - 2 v)} {(b + r) (b - r - 2 v)}) + 2 i \pi c_1}{\log(b - r) - \log(b + r)}$$

Using this I get a complex number for $$n$$ with real and imaginary part $$\not\in \mathbb{N}$$. If I compute $$x_n=q_n/q_{n-1}$$ and with this $$v$$ in $$\mathbb{R}$$ for an example it does work.
Any idea how to transform this to $$\mod P$$ (for values $$a,b$$ which don't have a root $$r$$)?
Or how to transform a $$v\in \mathbb{F}_P$$ to $$\mathbb{R}$$?

Update 4: Get rid of $$r$$ in $$q_n$$ equation
$$q_n$$ can be multiplied at top and bottom with $$r$$ to get $$r^2$$ only. With $$x0=0$$ the equation for $$x_n$$ would be (for $$n>2$$): $$x_n=q_n q_{n-1}^{-1}= \frac{2 \cdot \sum_{k=1}^{\lfloor n/2\rfloor} {n-1\choose 2k-1} b^{n-2k}(r^2)^{k} }{4\cdot \sum_{k=1}^{\lfloor (n-1)/2 \rfloor} {n-2 \choose 2k- 1} b^{n-2k-1} (r^2)^{k } } \mod P$$

The lower part need to be the inverse and not a divisor for $$\mod P$$.
However I don't see any way to extract a formula for $$n$$ out of this. Do I miss something?

• Have you actually tried evaluating your equation there mod $P$ with the finite-field variants of the operations (e.g. make sure to use multiplicative inverses rather than integer division) and found that it did not work? – Ella Rose Nov 17 at 19:49
• @EllaRose do you mean the equation for $q_n$? It does work if modulo used at the computed value $q_n$ before inserting it in the equation for $x_n$. Just using the multiplicative inverse don't work because the squareroot is $\in \mathbb{R}$. I thought it wont work but ty for reminding. The root can be calculated in $\mathbb{F}_P$ as well. Will check it out. – J. Doe Nov 17 at 21:00
• @EllaRose edited the question and added an equation mod $P$. However not for all values exists a root in $\mathbb{F}_P$. So far could not find any $a,b,P$ which produces all values in $\mathbb{F}_P$ and has a root for $4a+b^2$ – J. Doe Nov 17 at 21:57
• Does your closed formula for $q_n$ not work if you pick a root $r$ (i.e. $r$ with $r^2 = 4a + b^2$) over the extension field $\mathbb{F}_{p^2}$? – mk. Nov 19 at 22:59
• @mk. (not sure if I understood you correct) If I compute $q_{n}=a\cdot q_{n-2}+b\cdot q_{n-1} \mod P^2$ instead of $P$ I get the same results for $x_n$ (with $\mod P$). But that only works if a root $r$ exists. If that's the case the generator won't generate all numbers which is needed. In tests so far root $r$ with $r^2 = 4a+b^2 \mod P^2$ don't exist as well if there is no root $r$ with $r^2 = 4a+b^2 \mod P$ . ( For $a,b$ which don't generate all numbers I can find a root $\mod P^2$ and $P$, $q_n$ works with both of those) – J. Doe Nov 20 at 3:06

(too long for a comment)

You can calculate $$x_n$$ directly by your formulas, if you know, how to calculate in a quadratic extension field of a prime field $$\mathbb{F}_p$$.

Here is the code in python:

# global parameters
p = 65537
t = (p+1)/2
r = 313
s = 997

# from https://stackoverflow.com/a/9758173/99978
# modular inverse based on extended Euclidean algorithm
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
# for generating sequence x_n
def succ(x):
return (r*modinv(x, p)+s) % p

# check that 4*r+s^2 does not have any root mod p
for i in range(t):
if (i*i) % p == (4*r+s*s) % p:
raise Exception('not irreducible')

# conjugate in field
def conj(x):
a, b = x
return [a, -b]
# subtraction in field
def minus(x, y):
a, b = x
c, d = y
return [(a-c)%p, (b-d)%p]
# multiplication of a+bX with c+dX in field F_p/[X^2+4*r+s^2]
def mult(x, y):
a, b = x
c, d = y
return [(a*c+b*d*(4*r+s*s)) % p, (a*d+b*c) % p]
# exponentiation in field
def exp(x, n):
y = [1, 0]
while n > 0:
if n % 2 == 1:
y = mult(y, x)
x = mult(x, x)
n = n/2
return y
# inverse in field (same as exp(x, p*p-2))
def inv(x):
a, b = mult(x, conj(x))
if b != 0:
raise Exception('error in inverse')
return mult([modinv(a, p), 0], conj(x))
def qn(n):
x, y = mult(mult(minus(exp([s, 1], n), exp([s, -1], n)), exp([t, 0], n)), inv([0, 1]))
return x
# your formula for x_n given q_n
def xn(x0, n):
return ((qn(n+1)*x0+r*qn(n))*modinv(qn(n)*x0+r*qn(n-1), p)) % p

# a short test
x = 97
for i in range(500):
x = succ(x)
if xn(97, i+1) != x:
raise Exception('wrong result')


Extending a prime field $$\mathbb{F}_p$$ ("prime fields" are the minimal fields, the finite ones are the fields of order $$p$$, the infinite one is $$\mathbb{Q}$$) by the square root of an element $$z\in\mathbb{F}_p$$ that does not have a root $$\bmod p$$ works the same way as extending the real numbers $$\mathbb{R}$$ to the complex numbers $$\mathbb{C}$$: Just invent an element $$X$$ whose square is $$z$$ (for $$\mathbb{C}$$ this is usually called $$i$$), and calculate with numbers of shape $$a+bX$$ using $$X^2 = z$$.

Formally correct this can be done by calculation in the polynomial ring $$\mathbb{F_p}[X]$$ modulo the polynomial $$X^2-z$$ or in code like you see it above.

You can slightly optimize the code by using the fact that conjugation (replacing $$X$$ by $$-X$$) is commutes with addition/multiplication/exponentiation in the field, but I doubt that you will be able to invert the function $$n\mapsto x_n$$.

• well written code. That tighten the understanding of the root in an extension field. So there don't need to exist a complex number which can solve this but it behaves like this (+-*conj(x))). I also wrote some code for sequence generation including egcd, modinv before. To compute $n -> x_n$ comping $q_n$ in $\mathbb{R}$ and $\mod P$ also works. For planned usage other direction $x_n -> n$ should be hard to compute (~brute force). As you wrote I'm pretty sure I won't be able to invert this function :) Asked if it is only a problem for me or a general. Ty for new insights.Will have a look again – J. Doe Nov 25 at 23:05
• The sum of $x^n$ and its conjugate is just twice its "real" part, so $q_n$ is essentially the "real" part of an $n$-th power, but I doubt this will help for $x_n$. As you write that calculating over $\mathbb{R}$ (I assume double precision floating point numbers) gives you the correct results, your $p$ cannot be much longer than 40 bits. For $p$ up to around 64 bit, you can just generate a lookup table for $x_n\mapsto n$ using the python script with step size $2^{32}$ for $n$ and then iterate the generator for a given $x$ until you hit the table. – j.p. Nov 26 at 5:54
• Oh dear, how could I forget about this (lookup table hit). I though I finally found a solution for the target application :(. Will mark it as an answer. If you like you can add some hint about the lookup table in your main post just to be a complete answer to the question (in case other people interested in) – J. Doe Nov 26 at 16:13