Is there an efficient algorithm that allows to obtain a uniform distribution of all possible $4$-bit permutations from a single keyless $4$-bit S-Box?

Question

Let $S$ denote a keyless permutation that operates on $4$-bit inputs and returns $4$-bit outputs (that is, $S$ is a $4$-bit S-Box).

In this question, $x_0$ denotes an arbitrary bitstring that ends with 0000 (four zero bits) such that the length $L$ of $x_0$ is a multiple of $4$ and there is no upper bound for $L$. Then $x_0 + i$ denotes a bitstring that represents the sum of $x_0$ and $i$, assuming that $x_0$ is interpreted as a single natural number and all ($L - 4$) leading bits of $x_0$ are not affected. For example,

000001010000 + 1 =  000001010001,
000001010000 + 14 = 000001011110,
000001010000 + 15 = 000001011111

Does there exist an efficient iterative construction $F(x)$ based on $S$ so that the length of $F(x)$ is $4$ and a tuple $$\{F(x_0), F(x_0 + 1), \ldots, F(x_0 + 14), F(x_0 + 15)\}$$ is a sequence of $16$ different $4$-bit elements for any $x_0$ (that is, each $x_0$ generates a permutation of $4$-bit values), and all of the $2^4! = 20922789888000$ possible permutations are equiprobable if $x_0$ is a sequence of random bits followed by four zero bits (assuming that the length of $x_0$ is $4n$, where $n$ is an arbitrarily large natural number greater than $1$)? If it is not possible for all permutations to be absolutely equiprobable in the mathematical sense, then at least the distribution of permutations must be as uniform as possible and look “pseudo-random”.

Answer 1

In Haskell the recursive function

f s x = s (g 16 (x `div` 16) (x `mod` 16))
  where g 1 _ = id
        g n x = t (n-1) (x `mod` n) . g (n-1) (x `div` n)
        t x y z = if x == z then y
                            else if y == z then x else z

runs for any function s that is a permutation on [0..15] through all permutations on [0..15] if you let x run through all numbers 0...20922789888000*16-1.

Of course, s doesn't matter much, so I'd set s = id, and look at the definition of g using t: t x y is simply exchanging the values x and y, keeping all other values the same.

Now g n x is a(n arbitrary!) permutation on [0..n-1] depending on x: x mod n determines which element to exchange with n-1 (all other elements are untouched) and x divided by n is used as parameter to for defining recursively g (n-1).

If you pick n uniformly from 0..20922789887999, the permutation g n will be a uniform random element of the symmetric group on [0..15].