Fix a finite set of integers $S$ and a prime number $p$. Let $(a\_1, a\_2, \dotsc, a\_n)$, $(b\_1, b\_2, b\_3, \dotsc, b\_n)$ be two sequences of integers where the numbers $a\_i$ and $b\_i$ are chosen uniformly from the set $S$ (In the case I care about S = {-2,-1,1,2 }). I would like to understand what is the probability that $$\sum\_{i\_1 < i\_2 < \dotsb < i\_k } a\_{i\_1}a\_{i\_2}\dots a\_{i\_k} = 0 \mod p$$ an the probability that $$\sum\_{i\_1 \leq j\_1 < i\_2 \leq j\_2 < \dotsb < i\_k \leq j\_k} a\_{i\_1}b\_{j\_1}\dots a\_{j\_k}b\_{j\_k} = 0 \mod p.$$


I am interested in the asymptotic probabilities as $n$ increases and $k$ is fixed (or $k$ grows much slower than $n$) and expect some type of equidistribution as $n$. One particular case is [Question about estimating random symmetric sums modulo p](https://mathoverflow.net/questions/402034/question-about-estimating-random-symmetric-sums-modulo-p), where one particular sum is reduced to the linear case. Can this trick be for the general case?


Any suggestion, trick or vaguely related information related to this is highly appreciated.


**P.S.** Can the second type of sum (the one involving $a\_i$, $b\_i$) be somehow related to a sum of the first type?, the ordering makes things much more tricky and I am wondering whether there is a procedure to reduce this sum to the first type. (Also the second sum even though not so natural is the one that I really need to understand).


Thanks!



Distribution of some sums modulo p
Using the [Newton identities](https://en.wikipedia.org/wiki/Newton%27s_identities), one can (in the high characteristic regime $p>k$) express the elementary symmetric polynomial $\sum\_{i\_1 < \dots < i\_k} a\_{i\_1} \dots a\_{i\_k}$ in terms of the moments $\sum\_{i=1}^k a\_i^j$ for $j=1,\dots,k$. The question then boils down to the equidistribution of $\sum\_{i=1}^k (a\_i^j)\_{j=1}^k$ in ${\mathbf F}\_p^k$. There is however an obstruction to equidistribution because the points $(a^j)\_{j=1}^k$ for $a=-2,-1,1,2$ do not span the entirety of ${\mathbf F}\_p^k$ once $k \geq 4$, due to the existence of non-trivial polynomials $P$ of degree $4$ or higher that vanish at $-2,-1,1,2$. For instance because $P(x) = (x+2)(x+1)(x-1)(x-2) = x^4 - 5 x^2 + 4$ vanishes at these points, the vector $(m\_j)\_{j=1}^k := \sum\_{i=1}^k (a\_i^j)\_{j=1}^k$ is surely constrained to the hyperplane $m\_4 - 5 m\_2 + 4k = 0$. This is going to create some distortions to the probability that $\sum\_{i\_1 < \dots < i\_k} a\_{i\_1} \dots a\_{i\_k}$ vanishes mod $p$ that can be explicitly calculated for each $p,k$, but the formula is going to be messy (these errors will be of lower order in the limit $p \to \infty$ holding $k$ fixed though, by the Lang-Weil estimates).


(To put it another way, one can use the Newton identities to express the elementary symmetric polynomial as some polynoimal combination of the frequencies $d\_{-2}, d\_{-1}, d\_1, d\_2$ that count how often the random sequence $a\_1,\dots,a\_k$ attains each of its permitted values $-2, -1, 1, 2$. To count the probability that this polynomial vanishes mod $p$, one has to count the points in some variety over ${\bf F}\_p$ which in general is a task for the Lang-Weil estimate.)


Your second sum seems to be expressible in terms of a symmetric polynomial in a suitable matrix algebra. If one had $i\_k < j\_k$ in place of the constraint $i\_k \leq j\_k$ then one just needs to take the order $2k$ elementary symmetric polynomial of the matrices $\begin{pmatrix} 0 & b\_i \\ a\_i & 0 \end{pmatrix}$ and extract the top left coefficient. With $i\_k \leq j\_k$ the situation is more complicated but I expect there is still some sort of matrix representation. However being non-abelian I doubt there is a reduction to the abelian elementary symmetric polynomial considered earlier, and given how complicated that formula already was I'm afraid it is not going to be fun to try to control this sum (except possibly in the asymptotic regime where $p$ is somewhat large and $k$ goes to infinity; actually the nonabelian case might be substantially more "mixing" than the abelian one and one could conceivably get better asymptotics by using some of the theory of expansion of Cayley graphs, but this looks like a lot of work...).


