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

Question

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?

Answer 1

(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))
# your formula for q_n
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$.