CGMO, 2022

Let the set $D = \{0, 1, \dots, p-1\}$, and let the polynomial $f(x) = (x-1)x(x+1) = x^3-x$. The congruence and congruence symbol “$\equiv$” in this question refer to congruence modulo $p$. For $k=0, 1, 2, 3$, let
$$B_k = \{b \in D \mid \text{there are exactly } k \text{ elements } a \in D \text{ such that } f(a) \equiv b\}.$$
Consider an element $b$ in $B_2$, i.e., there exist $a_1 \not\equiv a_2$ such that $f(a_1) \equiv f(a_2) \equiv b$. In this case,
$$0 \equiv (a_1^3 - a_1) - (a_2^3 - a_2) = (a_1 - a_2)(a_1^2 + a_1a_2 + a_2^2 - 1) \Rightarrow a_1^2 + a_1a_2 + a_2^2 \equiv 1.$$
Thus, for $a_3 = -a_1 - a_2$, we have $a_3^2 + a_3a_1 + a_1^2 = (a_1 + a_2)^2 - (a_1 + a_2)a_1 + a_1^2 \equiv 1$, so it can be deduced that $f(a_3) \equiv f(a_1) \equiv b$, by the definition of $B_2$, $a_3$ is congruent to one of $a_1, a_2$, so either $a_2 \equiv -2a_1$ or $a_1 \equiv -2a_2$, assume the former, in this case, $a_1^2 + a_1a_2 + a_2^2 \equiv 3a_1^2 \equiv 1$, $3b \equiv 3a_1^3 - 3a_1 \equiv -2a_1$. When 3 is a quadratic residue modulo $p$, $a_1$ has two solutions; when 3 is not a quadratic residue modulo $p$, $a_1$ has no solution. Therefore, there are exactly $1 + (\frac{3}{p})$ solutions for $a_1$, and the corresponding $b$ also has $1 + (\frac{3}{p})$ solutions, i.e., $|B_2| = 1 + (\frac{3}{p})$.
$$\text{Obviously, } |B_1| + 2|B_2| + 3|B_3| = |D| = p.$$
On the other hand, consider the pairs $(u, v)$ that satisfy "$f(u) \equiv f(v)$ and $u \not\equiv v$" (called a collision). This is equivalent to $u^2 + uv + v^2 \equiv 1$ and $u - v \not\equiv 0$. After changing variables to $(x, y) \equiv (\frac{u-v}{2}, \frac{u+v}{2})$ (i.e., $(u, v) \equiv (x+y, y-x)$, such that $(x, y)$ and $(u, v)$ correspond one-to-one), it is transformed into $x^2 + 3y^2 \equiv 1$ and $x \not\equiv 0$.
We consider the number of $(x, y)$ pairs satisfying $x^2+3y^2 \equiv 1 \pmod p$ (i.e., $(3y)^2 \equiv 3-3x^2$).
$$\begin{aligned} T &= \sum_{x=1}^{p} 1 + \left(\frac{3-3x^2}{p}\right) = p + \left(\frac{-3}{p}\right) \sum_{x=1}^{p} \left(\frac{x^2-1}{p}\right) = p + \left(\frac{-3}{p}\right) \sum_{x=2}^{p} \left(\frac{\frac{x+1}{x-1}}{p}\right) \\ &= p + \left(\frac{-3}{p}\right) \sum_{x=1}^{p-1} \left(\frac{1+\frac{2}{x-1}}{p}\right) = p + \left(\frac{-3}{p}\right) \sum_{z=1}^{p-1} \left(\frac{z+1}{p}\right) = p - \left(\frac{-3}{p}\right). \end{aligned}$$
The above formula uses the fact that $z = \frac{2}{x-1}$ traverses the complete set of residues modulo $p$, and that $\sum_{k=1}^{p}(\frac{k}{p}) = 0$. Among the $T$ pairs $(x, y)$, there are exactly $1 + (\frac{3}{p})$ pairs for which $x \equiv 0$ (i.e., $3y^2 \equiv 1$). Therefore, the number of ordered collision pairs $(u, v)$ is
$$M = T - 1 - \left(\frac{3}{p}\right) = p - \left(\frac{-3}{p}\right) - 1 - \left(\frac{3}{p}\right) = T - |B_2|.$$
Thus, there are exactly $\frac{M}{2}$ unordered collision pairs $(u, v)$. Since the elements in $B_3$ correspond to 3 collisions, the elements in $B_2$ correspond to 1 collision, and the elements in $B_1$ and $B_0$ correspond to 0 collisions, it follows that $|B_2| + 3|B_3| = \frac{M}{2}$.
Thus, $|B_3| = \frac{M-2|B_2|}{6}$, so the number of all possible residues of $f(x)$ modulo $p$ is
$$|B_1| + |B_2| + |B_3| = p - |B_2| - 2|B_3| = p - \frac{M + |B_2|}{3} = p - \frac{T}{3} = \frac{2p + (\frac{-3}{p})}{3} = \left\lfloor \frac{2p + 1}{3} \right\rfloor.$$

Note 1: Another way to calculate the number of solutions $T$ is:
Consider the congruence equation $x^2 + 3y^2 \equiv z^2$, which is equivalent to $3y^2 \equiv z^2 - x^2 = (z+x)(z-x)$, having $p^2$ solutions ($y \equiv 0$ has exactly $2p-1$ solutions, $y \equiv 1$ has exactly $p-1$ solutions, and so on). The number of solutions with $z \equiv 0$ is
$$m_0 = 1 + (p-1)[1 + \left(\frac{-3}{p}\right)],$$
and the number of solutions for the remaining $z \equiv 1, z \equiv 2, ..., z \equiv p-1$ are equal (using $(x, y, z) \leftrightarrow (kx, ky, kz)$ to map the solutions with $z \equiv 1$ to the solutions with $z \equiv k$), each equal to $T = \frac{p^2-m_0}{p-1} = p - \left(\frac{-3}{p}\right)$.
If one is not familiar with Legendre symbols, one can first obtain $0 \le m_0 \le 2p$, and deduce that the number of solutions to $x^2 + 3y^2 \equiv 1$ is $T \in p-1, p, p+1$. The subsequent answer to this problem is $p - \frac{T}{3}$, which is an integer. From this, we can still get the needed answer.
Note 2: Another interpretation of the number of solutions $T = p - \left(\frac{-3}{p}\right)$ is: For each $y$, the number of $x$ that satisfy $x^2 \equiv 1 - 3y^2$ is $1 + \left(\frac{1-3y^2}{p}\right)$, therefore
$$\begin{aligned} T - p &= \sum_{y=0}^{p-1} \left( \frac{1 - 3y^2}{p} \right) \equiv \sum_{y=0}^{p-1} (1 - 3y^2)^{\frac{p-1}{2}} \\ &= \sum_{y=0}^{p-1} \left( \left(-3\right)^{\frac{p-1}{2}} y^{p-1} + \cdots + 1 \right) \equiv (-3)^{\frac{p-1}{2}} (p-1) = - \left( \frac{-3}{p} \right). \end{aligned}$$