Find all the primes $p$ so that $$x^3+x \equiv 1 \pmod p \tag{1}$$ has integer solutions.
We consider $x$ and $y$ are the same solutions iff $x\equiv y \pmod p.$
We can prove that $(1)$ cannot have exactly $2$ solutions unless $p=31$. When $p=31$, the solutions of $(1)$ are $x \equiv {17\text{ (double), } 28} \pmod {31}.$
My question is: when does $(1)$ have no solution, and when does $(1)$ have $1$ solution and when does $(1)$ have $3$ solutions? (As Ma Ming has pointed out, an equation of degree $3$ has no more than $3$ solutions.)
I have proved that if $p \not = 31$, and $a$ is a solution of $(1)$, then $(1)$ has $3$ solutions iff $$\left(\frac{a-1}{p} \right)=\left(\frac{a+3}{p} \right),$$ here $\left( \frac{a}{p} \right)$ is the Jacobi symbol. Thanks in advance!
To augment the answers of Pete Clark and Will Jagy: the splitting field for $x^3 + x - 1$ is the Hilbert Class Field of $\mathbb Q(\sqrt{-31})$. As Pete noted, its Galois group over $\mathbb Q$ is isomorphic to $S_3$. Thus the splitting behaviour of a prime $p$ depends first on whether is is a square or not mod $31$, and (if it is a square) whether or not it splits principally in $\mathbb Q(\sqrt{-31})$; this is the theoretical origin of the quadratic forms appearing in Will Jagy's answer.
If we embed $S_3$ into $GL(2,\mathbb C)$ in the usual way, then we will get a two-dimensional Galois representation whose splitting field is precisely the splitting field we are discussing, and it will correspond (via the modern interpretation of results of Hecke) to a modular form on $\Gamma_1(31)$ of quadratic nebentypus, which one can write down as a difference of theta functions (associated to the quad. forms in Will Jagy's answer).
So there is a certain $q$-expansion $\sum_{n = 1}^{\infty} a_n q^n$ with integer coefficients $a_n$ such that $x^3 + x - 1$ splits mod $p$ (different from $31$) if and only if $a_p = 2$. (A priori, $a_p = 0,-1$, or $2$.)
If you had asked the corresponding question for $x^3 - x + 1$, whose discriminant equals $-23$, then the same story would apply (with $31$ replaced by $23$ everywhere), and I could have given you the following formula for the $q$-expansion, namely $q \prod_{n = 1}^{\infty} (1-q^n)(1-q^{23 n})$. E.g. in this case the first $p$ such that $a_p = 2$ (and hence such that $x^3 - x + 1$ splits mod $p$) is $p = 59$.
In the case of $-31$, though, I don't know a simple formula for the $q$-expansion, although you can certainly compute any number of terms of it that you want using modular forms software.