Solve $x^3+x \equiv 1 \pmod p$ 
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!
 A: The discriminant of the irreducible integer polynomial $f(x) = x^3+x-1$ is $-31$.  Thus (as you've said) there will be repeated roots iff $p = 31$.
Let $F = \mathbb{Q}[x]/(x^3+x-1)$.  The cubic number field $F$ is not Galois; thus its Galois closure $M$ is an $S_3$-extension of $\mathbb{Q}$.  By basic algebraic number theory, the splitting pattern of $f(x)$ modulo $p$ corresponds to the splitting pattern of the prime ideal $(p)$ in $F$.  For instance we have three roots iff $p$ splits completely, and we have exactly one root iff $p \mathbb{Z}_F = \mathfrak{p}_1 \mathfrak{p}_2$, where $\mathfrak{p}_1$ is a degree one prime and $\mathfrak{p}_2$ is a degree two prime.  
Further, by the Cebotarev Density Theorem one can compute the densities of these sets of primes by passing to the Galois closure $M$.  A prime splits completely in a number field iff it splits completely in its Galois closure -- equivalently every Frobenius element over $p$ is trivial, and by Cebotarev Density the set of such primes has density $\frac{1}{6}$. A prime $p \neq 31$ which remains inert in $F$ must split into two degree $3$ primes in $M$, since otherwise we would have an order $6$ Frobenius element over $p$ and $S_3$ contains no element of order $6$.  By Cebotarev, the density of the set of such primes is equal to the proportion of $3$-cycles in $S_3$, so $\frac{1}{3}$.  The last case is (just) a little messier than the others, but we can skip it since the densities must add to $1$.  Thus:
$\bullet$ If $p = 31$, $f$ has a multiple root modulo $p$.  Otherwise it has distinct roots in $\overline{\mathbb{F}}_p$.  
$\bullet$ The set of primes $p$ such that $f$ has three roots modulo $p$ has density $\frac{1}{6}$.  
[Note that the more naive guess would have been $\frac{1}{3}$.  One really does have to pass to the Galois closure to apply Cebotarev!]
$\bullet$ The set of primes $p$ such that $f$ has one root modulo $p$ has density $\frac{1}{2}$.  
$\bullet$ The set of primes $p$ such that $f$ has no roots modulo $p$ has density $\frac{1}{3}$.  
Now you may ask why I am only telling you the densities of these various sets rather than telling you more explicitly which primes fall into which.  The answer is that when the Galois group of the polynomial is nonabelian -- as here -- there will in general be no simpler description of the primes than the above Cebotarev conditions.  In contrast, when the polynomial has abelian Galois group the sets of primes can simply be given by congruence conditions.  For instance, I randomly found this handout which treats the similar looking polynomial $x^3-3x-1$ instead, but this polynomial has cyclic Galois group, which makes all the difference in the world.
It is true that in certain circumstances one can find explicit descriptions on the primes which are a bit more complicated than just congruence conditions.  Sometimes they come out in terms of higher residues and sometimes they come out in terms of modular forms!  But to the best of my knowledge (and I wish I had more knowledge here) this is rather sporadic and lucky; off the top of my head I'm not sure that this kind of nice description exists here.
A: 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.  
