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
$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?

  • $\begingroup$ 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? $\endgroup$
    – Ella Rose
    Nov 17, 2019 at 19:49
  • $\begingroup$ @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. $\endgroup$
    – J. Doe
    Nov 17, 2019 at 21:00
  • $\begingroup$ @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$ $\endgroup$
    – J. Doe
    Nov 17, 2019 at 21:57
  • $\begingroup$ 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}$? $\endgroup$
    – mk.
    Nov 19, 2019 at 22:59
  • $\begingroup$ @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) $\endgroup$
    – J. Doe
    Nov 20, 2019 at 3:06

1 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)
    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')
    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$.

  • $\begingroup$ 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 $\endgroup$
    – J. Doe
    Nov 25, 2019 at 23:05
  • 1
    $\begingroup$ 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. $\endgroup$
    – j.p.
    Nov 26, 2019 at 5:54
  • $\begingroup$ 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) $\endgroup$
    – J. Doe
    Nov 26, 2019 at 16:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.