I'm looking for functions $f_1, f_2$ which are commutative to each other.
$$f_1(f_2(x)) = f_2(f_1(x))$$
Random numbers $r_{i,j}$ should be computed in a chain/sequence:
$$r_{i+1,j} = f_1(r(i,j)) $$
$$r_{i,j+1} = f_2(r(i,j)) $$
$$r_{i+1,j+1} = f_1(f_2(r(i,j))) = f_2(f_1(r(i,j)))$$
They need to produce many but finite amount of numbers. If with given $r_{0,0}$ the funtion $f_1(x)$ can produce $s_1$ unique numbers and $f_2(x)$ can produce $s_2$ unique numbers than $f_1(f_2(x))$ should be able to produce $N$ numbers with $N>s_1,s_2, s_1+s_2$. It should be at least $\mathcal{O}(s_1\cdot s_2)$
There should be no close form for $r_{i,j}$ with given $r_{0,0}$. Or as alternative the computation as least as hard as generating all variables in between.
If a random variable $v$ is known to be part of the sequence starting at $r_{0,0}$ with $v = r_{a,b}$ for unknown $a,b$ it should take in mean $\mathcal{O}(a\cdot b \cdot \mathcal{O}(f) )$ operations to find $a,b$. That also means if $a,b$ is known it does not help for finding another sequence element $r_{c,d}$ for the case $c=a$ or $d=b$.
EDIT: As Misha Lavrov wrote in comments $\mathcal{O}(a\cdot b \cdot\mathcal {O}(f) )$ is not possible if there is also a inverse function with same computation time. Only $\mathcal{O}(a + b)\mathcal{O}(f)$ steps are needed. For simplicity I only wrote about two functions $f_1,f_2$. In target usage I have 3 of them so it would be $r_{a,b,c}$. For this the product should work again at least with two factors. With 3 functions $\mathcal{O}( (a + b)\cdot c)\mathcal{O}(f)$ steps are needed ($a,b,c$ can be permuted to get lowest value).
(end of EDIT)
If $f_1,f_2$ is not valid for all $r_{0,0}$ there should be a way to validate if a random value $r_{0,0}$ is suitable or not.
Furthermore there need to exists a inverse function of $f_1, f_2$ which has about the same computation time.
Solution can include some special kind of groups, rings, hypercomplex numbers or anything else. Only condition $r_{i+1,j}, r_{i,j+1}$ need to be computable at a PC with storage less than the total number of unique numbers $N$.
Target application
Needed for an arrangement of unique values with certain size. Too many to store them all, too less to use common cryptography methods like elliptic curves in each dimension. Only the product of the size in each dimension is big enough to have some long enough computation time. Use case need 3 functions $f_1,f_2,f_3$ which are commutative. The randomness of those values don't need to pass any tests for randomization. The only requirement is they give no clue about the index related to $r_{0,0}$.
The $r_{0,0}$ is random. Same $r_{0,0}$ should give same set $R=\{r_{i,j} \forall i,j\}$. A different $r_{0,0}'$ can produce a different set $R'=\{r_{i,j}' \forall i,j\}\neq R$. The amount of different sets is limited. It should be much smaller than the amount of different $r_{0,0}$
Functions which don't work
-Binary XOR functions are commutative but they can only produce a very small amount of unique numbers.
-Cyclic group with generators $g,h$ and suitable prime $P$ which has the form:
$$g^i h^j \mod P$$
$$f_1(x) = x\cdot g \mod P$$
$$f_2(x) = x\cdot h \mod P$$
This can be reduced to a problem of just finding $i$ or $j$. So the problem scales with the max cycle size and not with their product.
-Chebyshev polynomials are commutative functions. Did some tests with $$f_1(x) = T_n(x) \mod P$$ $$f_2(x) = T_{n+1}(x) \mod P$$ with suitable $P$ and $r_{0,0}$ they can produce up to a $P/4$ different numbers. But if $f_1,f_2$ used together the total amount of different numbers won't increase. For same $r_{0,0}$ they produce same numbers with just a different order. (They also need some settle in time but that's ok)