How did we discover that this quadratic residue oriented PRNG generates unique numbers in a sequence?

Question

Question (tl;dr)

How do we know that even for extremely large numbers like even far past $32^{15}$ (any bigint, or any number really), that as you increment the sequence from $0$ to $n$, $n$ being the max, that every output as a result of permute(i) (defined below), for every i up to n, will be unique? How do we know that? Who discovered it? (not necessary for the answer to know exactly who discovered it, but it would be a bonus to know!)

To me the only way to show this would be to check the output of permute(i) has no duplicates running through every possible input (up to n). But how do we know how it works without doing that? How did anyone ever figure this out? (When, roughly speaking, was this even discovered)?

Background

There is a "quadratic residue oriented" PRNG "magic" function (to me it's magic), described here, an example of which is found below, which I got from asking this SO question and made into a JavaScript library. It can be implemented in just 2 lines of code (or 10):

const rotate = (x, o) => x >= o ? x : (x <= (o / 2n)) ? ((x * x) % o) : o - ((x * x) % o)
const permute = (v, w, s, x, m) => rotate((rotate(v, w) + s) % m ^ x, w)

Sidenote: I don't know if I am naming these functions correctly, I just called them what seemed appropriate.

An example of how it is used is like this:

const E = 3132343537383103113163n
const A = 975319753197531975319n
const O = 541613713n
const U = 32 ** 15

let i = 0
while (i < U) {
  const x = permute(i, E, A, O, U)
  console.log(x)
  i++
}

More description is in the repo description, on why we chose specific primes and a use case for this function.

Now, the only information I have on how this works comes from that blog linked above, which is copied here:

---

---

Conclusion

How is this even possible? Seemingly random numbers jumping around an enormous sequence, fill in the slots without ever encountering a duplicate?! That is so crazy to me. Hence the main question.

Answer 1

As the linked blog explains, as long as each operation you do just shuffles the integers around, you'll never get duplicates. You can do as many shuffles of a pack of cards as you like, you'll never get duplicate cards appear when you deal them out later.

Shuffle operations on a finite set are also called permutations. They are invertible - any shuffle could be undone by the appropriate unshuffling operation.

The three shuffling operations used here are:

• addition of $s$ modulo $m$
• xor-ing with x
• the operation which you call "rotate"

The first two are permutations of $0$...$m-1$ because they are easily invertible (assuming $m$ is a power of $2$). It is much trickier to see why the third operation is a permutation. Note that this third operation works if $w$ is a prime which satisfies $w\equiv 3 \pmod 4$ but not necessarily in other cases.

There are a few steps needed to prove this.

1. The non-zero squares mod $p$ cover exactly half of the numbers $1$,...,$p-1$.

Proof: The numbers $1^2$,...,$(p-1)^2$ cover at most half the numbers cause every number is duplicated due to the fact that $x^2=(-x)^2$. However, there can't be further duplicates:

$$x^2\equiv y^2 \pmod p \\ x^2-y^2\equiv 0 \pmod p \\ (x-y)(x+y)\equiv 0 \pmod p \\ x-y\equiv 0 \pmod p \text{ or } x+y\equiv 0 \pmod p \\ x\equiv y \pmod p \text{ or } x\equiv -y \pmod p$$

Note that we used the fact that $p$ is a prime in the step to the second to last line, using the fact that if a prime divides a product, then the prime must divide one of the factors.

2. If $-1$ is not a square mod $p$, then the squares and the negated squares cover all of the numbers $1$,...,$p-1$.

Proof: A non-zero square cannot equal a negated square, because if $x^2\equiv -y^2 \pmod p$ then $(xy^{-1})^2 \equiv -1 \pmod p$, making $-1$ a square contrary to our assumption. This shows that the non-zero squares and negated squares don't overlap. They each cover half the numbers, so together they cover all the numbers.

Putting this together, you see that the squares of $1$,...$(p-1)/2$ and the negated squares of the rest cover all the non-zero numbers mod $p$. Combined with mapping $0$ to $0$, you get a permutation of the first $p$ integers. We do need $p$ to be a prime such that $-1$ is not a square modulo $p$. This is the case when $p\equiv 3 \pmod 4$, which can be shown using Euler's criterion.

Addendum:
Although this "rotate" operation is a permutation of the numbers $0$ to $p-1$, and therefore invertible, it is not easy to invert when $p$ is large. Suppose the result of the operation is $r$, and the input we want to find is $x$. It is relatively easy to check whether $r$ is a square mod $p$ using quadratic reciprocity. If $r$ is a square, then $x=\sqrt{r} \pmod p$, and otherwise $x=-\sqrt{-r} \pmod p$. The hard part is taking that square root mod $p$, but there are algorithms for that such as the Tonelli-Shanks algorithm.