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, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.