Question about estimating random symmetric sums modulo p Let $n > 0$ be a positive integer (large) and $p > 2$ a fixed prime number. What is the probability that $$\sum_{ 1 \leq i < j \leq n} a_ia_j = 0 \mod p$$ where $a_1, a_2, \dots a_n$ are chosen uniformly from the set $S = \{-1, 1\}$. Does this sum equidistribute mod $p$ as $n$ goes to infinity? What would be the speed of equidistribution in terms of $n$? Is there any literature in this type of random sums? I would be surprised if not but I am unable to find anything related or similar to this.
One can also ask what is the probability of this sum being actually zero, but I also have no idea how to deal with it and thought that modulo a prime would be simpler.
 A: The condition
$$\sum_{ 1 \leq i < j \leq n} a_ia_j \equiv 0 \pmod p$$
is equivalent to
$$\left(\sum_{ 1 \leq i\leq n} a_i\right)^2 \equiv n \pmod p.$$
So a necessary condition is that $n$ is a quadratic residue modulo $p$ (including the zero residue). If $n$ is divisible by $p$, then the above condition says that the sum of the $a_i$'s is divisible by $p$. Otherwise, the condition says that the sum of the $a_i$'s is congruent to one of the two square-roots of $n$ modulo $p$. Now it is easy to see that the sum of the $a_i$'s is equidistributed modulo $p$ (think about what happens when an $a_i=1$ is switched to $a_i=-1$), hence in the first case the probability is $1/p+o(1)$, in the second case it is $2/p+o(1)$, as $n$ tends to infinity.
In fact the  probabilities can be calculated explicitly as a linear combination of $n$-th powers of $p$ complex numbers (which only depend on $p$), since the sum of the $a_i$'s modulo $p$ is determined by $\#\{i:a_i=1\}$ modulo $p$, and vice versa. Compare with this post, where the role of $p$ is played by $4$. It follows, in particular, that the $o(1)$ terms above decay exponentially fast. For a more complete reference, see Theorems 8.7.2 & 8.7.3 in Wagner: A first course in enumerative combinatorics (AMS, 2020).
