Any speed up methods for finding the index of a random value produced by the Inversive congruential generator?
(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 &gt; 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$.