I find this problem intriguing. The following is based on the ideas in the comments, so please join me in thanking Aphelli.
As per Igor's initial work the polynomial $f(x)$ has discriminant $-18^2$. It is well known that an irreducible cubic over a finite field has a cyclic Galois group, hence its discriminant must be a square. This implies that $-1$ is a quadratic residue modulo $p$, so $p\equiv1\pmod4$. We will be needing information about $p$ modulo three, so we split the handling of this problem into two cases: $p\equiv 1\pmod{12}$ and $p\equiv 5\pmod{12}$.
Cardano's formula tells us that the roots $\alpha,\beta,\gamma$ have the form
$u+v$, $\omega u+\omega^2v$, $\omega^2u+\omega v$, where $u$ and $v$ satisfy
$uv=-(-3/3)=1$ and $u^3+v^3=-4$. That system has solutions
$$
u^3=-2+\sqrt3,\qquad v^3=-2-\sqrt3.\tag{1}
$$
Here $\omega$ is a primitive third root of unity, and $\sqrt3$ obviously is a root of the equation $x^2=3$. A priori both of those reside in some extension field of $\Bbb{F}_p$ as do $u$ and $v=1/u$.
Assuming first that $p\equiv1\pmod{12}$.
In this case quadratic reciprocity tells us that there exists a square root of three in the field $\Bbb{F}_p$. Assume that we have located one such. Let's fix it once and for all. Let us select an element $u$ such that $u^3=-2+\sqrt3$. Rewriting Cardano's formula tells us now that
$$
\begin{aligned}
\alpha&=u+\frac1u,\\
\beta&=u\omega+\frac1{u\omega},\\
\gamma&=u\omega^2+\frac1{u\omega^2}.
\end{aligned}
$$
In the irreducible case the roots are not in the prime field, so they are each others images under Frobenius. Indeed, we can then dictate $\alpha^p=\beta$, $\beta^p=\gamma$, because we have not yet specified $\omega$! However, we need a way to figure out whether these elements are in the prime field, and also a way to select the correct $\omega$.
Because $-2+\sqrt3$ is in the prime field, we have $(-2+\sqrt3)^p=-2+\sqrt3$. In other words, we have the equation
$$
u^{3p}=(-2+\sqrt3)^p=-2+\sqrt3=u^3.
$$
So if we define
$$z:=u^{p-1}=(-2+\sqrt3)^{(p-1)/3},\tag{2}$$ it follows that $z^3=1$.
Observe that it is possible to have $u^p=u$, or equivalently $z=1$.
In this case $u$ is an element of the prime field, and so is the root $\alpha$, taking us back to the non-interesting case. In the interesting cases $z\neq1$, and we can declare
$$
\omega=z=(-2+\sqrt3)^{(p-1)/3}
$$
to be the key primitive third root of unity. This has the important consequence that
$$u^p=\omega u,\tag{3}$$
implying that $\alpha^p=\beta$, $\beta^p=\gamma$ and $\gamma^p=\alpha$.
So the upshot is that with this particular choice of $\omega$ the Frobenius then acts as $\alpha\mapsto\beta\mapsto\gamma\mapsto\alpha$.
We, finally, arrive (after a straightforward calculation assisted by Mathematica) at the formula
$$
P:=(\beta-\alpha)(\gamma-\beta)(\alpha-\gamma)=3\omega(1-\omega)\frac{u^6-1}{u^3}.
$$
Given $u^3=-2+\sqrt3$ it follows that $(u^6-1)/u^3=-2\sqrt3$, so
$$
P=-6\sqrt3\omega(1-\omega). \tag{4}
$$
A point here is that as long as we picked $\omega$ according to $(3)$, then the familiar equation
$$\omega-\omega^2=i\sqrt{3}\tag{5}$$
gives the proper sign for $i$! The choice of $\sqrt3$ is immaterial, because only $(\sqrt3)^2$ appears, when $i$ is selected according to $(5)$.
An example case of the possibility $z=1$ is $p=97$. Then we can choose $\sqrt3=10$, $-2+\sqrt3=8$ and $8^{(p-1)/3}=2^{p-1}\equiv1\pmod p$. So in this case $f(x)$ has roots in the prime field (a case Igor solved himself). It also follows that in this case $f(x)$ splits completely over the prime field, because the other roots are then also fixed points of Frobenius. It would be nice to know for which primes $p$ this happens. I believe class field theory has something to say about this problem, but I don't have the time to dig into that. Anyway, calculating $z=(-2+\sqrt3)^{(p-1)/3}$ is a tool to diagnose this possibility.
All of the above can be utilized in many ways to tackle the problem. Here's how I would go about it. I make no claims about this being the cleanest way, but it has a reasonable computational complexity (assuming you are familiar with basic tricks of the trade).
Assume $p\equiv1\pmod{12}$.
- Pick random elements $a$ of $\Bbb{F}_p$. Calculate $\tilde{\omega}=a^{(p-1)/3}$. If $\tilde{\omega}=1$, discard it, and try another choice of $a$. After this step we have a "proto-$\omega$".
- Similarly find a "proto-$i$", denote it $\tilde{i}$, by raising random elements to power $(p-1)/4$ until you get something $\neq\pm1$.
- The point of steps 1 and 2 was to give a fast way of locating a modular square root of three. This comes from the relation (of complex numbers)
$$\omega=\frac{-1+i\sqrt3}2\implies \sqrt3=\frac{2\omega+1}{i}.$$
Use the proto-versions $\tilde{\omega}$ and $\tilde{i}$ on the right hand side to arrive at an element of the prime field that can take over the duties of $\sqrt3$. It may well be that computational number theory knows of a more efficient way of finding a modular square root of $3$. If so, then the preliminary steps 1 and 2 can be left out. Anyway, at this point we can fix a choice of $\sqrt3\in\Bbb{F}_p$.
- Calculate
$$z=(-2+\sqrt3)^{(p-1)/3}.$$
If $z=1$ conclude that the answer is equal to zero (as $f(x)$ factors). Otherwise, define $\omega:=z$.
- Calculate
$$P:=(\alpha^p-\alpha)(\beta^p-\beta)(\gamma^p-\gamma)=-6\sqrt{3}\omega(1-\omega).$$
An account of how the case $p\equiv5\pmod{12}$ differs from the former case.
Cardano's formula still works as above. We also have $\sqrt{-1}\in\Bbb{F}_p$, but $\sqrt3$ and consequently also $\omega$ are in the quadratic extension field $\Bbb{F}_{p^2}$. A consequence of this is that the Frobenius will map $u^3=-2+\sqrt3$ to its conjugate $-2-\sqrt3=v^3=1/u^3$.
Therefore this time
$$u^p=1/(uz)\tag{6}$$
for some $z$ that satisfies $z^3=1$. Continuing the analogy, this means that $\alpha^p=\beta$ assuming that $z\neq1$ and that we select
$$
\omega=z=\frac1{u^{p+1}}=\frac1{(-2+\sqrt3)^{(p+1)/3}}.
$$
Observe how nicely the congruence arithmetic fits into this picture as this time $p+1$ will be divisible by three! This is, of course, not a coincidence but rather a well understood property of elementary Galois theory as well as the structure of finite fields. This time only $p^2-1$ is divisible by three, so the third roots of unity take us to the quadratic extension.
The expansion for
$$
P=(\beta-\alpha)(\gamma-\beta)(\alpha-\gamma)=-6\sqrt3\omega(1-\omega)
$$
should (?) be the same as above (IIRC the derivation of formula $(4)$ only needed $\omega^3=1$ and $u^3=-2+\sqrt3$). Again, if $z=1$, then $\alpha^p=\alpha$ and the polynomial $f(x)$ will be reducible.
So this case turned out to be very similar. However the above calculations take place in the field $\Bbb{F}_p[\sqrt3]$, so to carry out the same algorithm you need to implement the arithmetic of that field as well.
