What are the steps involved in solving a quartic polynomial modulo a prime modulus? This:
$$x^4 + 21x^3 + 5x^2 + 7x + 1 \equiv 0 \mod 23$$
Leads to:
$$x = 18 || x =19$$
I know this because of this WolframAlpha example and because a fellow member posted it in a since deleted & related question. 

What I don't understand are the steps involved in arriving at x = 18 || x = 19 from this equation.
My question starts with the reduced terms mod 23 example in the linked question. I'm now trying understand how to reduce this equation to x = 18 || x = 19.
I have come across a few posts and theorems that hint a solutions, but I lack the math skills to connect any of it together. I am a software developer, not a mathematician. So if anyone can walk me through some steps on how to get from the equation to 18 || 19, that would be great! 
This is a toy example representing a new Elliptic Curve Crypto operation where the actual modulus is $2^{256}$ large. So, trying all possible values x is not practical. WolframAlpha is capable of producing solutions to my large modulo equations in a fraction of a second so I know they aren't trying all possible values x.
Fermat’s Little Theorem seems the most promising so far, but I don't understand how to apply it to this equation. This post describes a solution but unfortunately their example is very basic and not very relatable to my equation.
Anything would be helpful here. Steps would be great. Thanks!
 A: If I were asked to "solve" a (monic, integer) quartic polynomial modulo a prime modulus ($23$ in the toy problem described here), I would first determine whether the polynomial can be factored over the rationals (equiv. over the integers by Gauss's lemma).

Here the polynomial turns out to be irreducible over the integers:
  $$ f(x) := x^4 + 21x^3 + 5x^2 + 7x + 1 $$

If there were a factor of degree one in $\mathbb Z[x]$, then by the Rational Roots Theorem there would be a root $\pm 1$.  It is easily checked that this is not the case.  The only other possible factorization over $\mathbb Z[x]$ would be the product of two quadratics:
$$ (x^2 + ax + 1)(x^2 + bx + 1) $$
or:
$$ (x^2 + ax - 1)(x^2 + bx - 1) $$
These possibilities can be ruled out by comparing the coefficients of $x^3$ and $x$ that would result, since this gives inconsistent values of $a+b$.
It is a minor frustration, but if $f(x)$ did factor over the integers, it would also factor over the integers mod $p=23$.  The converse is not valid.  It often happens that polynomials factor modulo an integer but are irreducible over the rationals (integers).

We now come to a connection with Fermat's Little Theorem:
  $$ x^p \equiv x \bmod p $$
  for any prime modulus $p$.  

Not only are all residues $a = 0,1,\ldots,p-1$ mod $p$ roots of $x^p - x$, this $p$th degree polynomial is the exactly the product of all $p$ of the first degree irreducible polynomials mod $p$.  See these class notes (Prop. 1) for a more general proposition for all finite fields.
We proceed to compute the polynomial GCD of $f(x)$ and $x^p - x$, which will give us the product of any first degree factors of $f(x)$.  If $f(x)$ splits over the integers mod $p$ (factors completely into first degree polynomials), we would get $\gcd(f(x),x^p-x)=f(x)$ back.  That would mean $f(x)$ has four distinct roots without telling us what those are!  But in the present case (with two distinct roots), we will instead get $f(x)$ factored as a product of two quadratics mod $p$.
Our chances of getting distinct factors are improved somewhat by noticing how easily factored $x^p - x$ is for odd primes $p$:
$$ x^p - x = x\left(x^{\frac{p-1}{2}} + 1\right)\left(x^{\frac{p-1}{2}} - 1\right) $$
Thus, instead of calculating $\gcd(f(x),x^p-x)$ we can calculate the GCD of $f(x)$ with each of those (coprime) factors of $x^p-x$.  This gives a chance of finding a first degree factor in one place and another first degree factor in another place.
By inspection we see that $\gcd(f(x),x) = 1$ because the constant term of $f(x)$ is nonzero.  Now with $p=23$ the two interesting factors of $x^p-x$ become $x^{11}+1$ and $x^{11}-1$.  We will compute both of their GCD's with $f(x)$, and as it turns out, we will get both of the two distinct first degree factors that way.
Since $x^{11}$ is a "shared" intermediate result, we compute its remainder modulo $f(x)$ and save the effort of doing that twice.  It turns out:
$$ x^{11} \equiv 9x^3 - 8x^2 - 2x + 5 \bmod{f(x)} $$
So the first step in finding $\gcd(f(x),x^{11}+1)$ is getting the remainder of $x^{11}+1 \bmod f(x)$ is $9x^3 - 8x^2 - 2x + 6$.  Note that we needed to preserve the nonmonic leading term of $x^{11} \bmod f(x)$ because we had to add $+1$ (resp. $-1$) correctly.
However for the follwing steps of the Euclidean algorithm for polynomials it is allowed to factor out that leading coefficient and work only with monic polynomials as the divisors:
$$ 9x^3 - 8x^2 - 2x + 6 \equiv 9(x^3 - 6x^2 + 10x - 7) \bmod 23 $$
Thus the next "division algorithm" step gives us:
$$ f(x) \equiv (x+4)(x^3 - 6x^2 + 10x - 7) - 4x^2 - 3x + 6 \bmod 23 $$
The remainder here becomes our divisor in the next step, so normalizing:
$$ -4x^2 - 3x + 6 \equiv -4(x^2 - 5x + 10) \bmod 23 $$
And so we continue the Euclidean algorithm:
$$ x^3 - 6x^2 + 10x - 7 \equiv (x-1)(x^2 - 5x + 10) - 5x + 3 \bmod 23 $$
$$ -5x + 3 \equiv -5(x+4) \bmod 23 $$
$$ x^2 - 5x + 10 \equiv (x-9)(x+4) + 0 \bmod 23 $$
This last remainder being zero tells us the GCD is found:
$$ \gcd(f(x),x^{11}+1) = x+4 $$
As a first degree factor of $f(x)$, this identifies one of its roots is $-4$ or equivalently modulo $23$, $x=19$.
A similar computation gives $\gcd(f(x),x^{11}-1) = x+5$, which identifies the other roots as $-5$ or $x=18 \bmod 23$.
Because $p=23$ was asked as a "toy problem", I'll point out two ways that computing with a large prime affects the complexity of factoring a quartic polynomial over that field of coefficients.  (to be continued)
A: The OP requested that I link my other answer as an answer to this one as well. 
A: I too believe, like saulspatz, that for small moduli one might just try all the possible values.
Another idea that might work for some simple equations is the following one, although it should be a last-resort technique (here I was able to make it work only because I already knew the solutions):
Since
$$21 = -2 + 23,\quad 5 = -64 + 3 \cdot 23, \quad 7 = -85 + 4 \cdot 23, \quad 1 = 300 - 13 \cdot 23$$
the equation is equivalent to:
$$x^4 - 2 x^3 - 64x^2 - 85 x + 300 \equiv 0 \pmod {23}$$
Now, by the integral root theorem, we check whether some divisors of $300$ are roots of the polynomial over $\mathbb Q$. Indeed,
$$(-4)^4 - 2(-4)^3 - 64 (-4)^2 - 85 (-4) + 300 = 0$$
$$(-5)^4 - 2(-5)^3 - 64 (-5)^2 - 85 (-5) + 300 = 0$$
We divide the polynomial by $(x + 4)$ and $(x + 5)$, obtaining:
$$(x + 4)(x + 5)(x^2 - 11x + 15) \equiv 0 \pmod {23}$$
Finally, since $\Delta = (-11)^2 - 4 \cdot 15 = 61 \equiv 15 \pmod {23}$ and $15$ is not a quadratic residue modulo $23$, the only solutions are $-4$ and $-5$.
A: If one of the roots ($x=19$) is known then the decomposition of the equation  isn't hard.
The substitution
$$x=y-4,\tag1$$
gives the least sum of the coefficients, wherein one of the roots must be zero:
$$y^4+5y^3-151y^2+719y-1035=0,$$
$$y^4+5y^3+10y^2+6y\equiv0\pmod{23}.\tag1$$
If the roots are not known then the easiest way is checking of polynomial values by modulo $23$.
The Vieta theorem can increase the bust by the next way.
If $x=0,$ then the polynomial value is 1, with the divisors $\pm1.$
If $x=1,$ then the polynomial value is 12, with the new divisors $\pm2, \pm3, \pm4, \pm6, \pm12$ etc.
This allows to check only the possible values.
The equation $(1)$ can be decomposed in the form of
$$y(y+1)(y^2+4y+6)\equiv0\pmod{23},\tag2$$
with the roots $y\equiv-1,0\pmod{23},$
$$\mathbf{\color{brown}{x\equiv18,19\pmod{23}.}}$$
The equation becames cubic. The previous way can be used.
At the same time, the quadratic equation
$$y^2+4y+6\equiv 0\pmod{23}$$
is well known. It does not have the integer roots.
In partial, this can be proved, using the tables of quadratic residues. But if the modulo is small, the bust seems easier.
A: The general methods for solving quartics via radicals work modulo $23$; IIRC they work in every characteristic except for 2 and 3, so you could apply those if you know how to take square and cube roots. This will often require some intermediate calculation in extension fields
23 is a small so simply trying every possible value and checking if it is a root is feasible, especially via program. Of course, this is less feasible for large primes.
The general method for this sort of problem, however, is basically to apply a general polynomial factorization algorithm for finite fields to discover the linear factors of your polynomial.
The fact you're just looking for the roots rather than the full factorization doesn't really simplify these general methods, although with care it would let you do less work. For example, if you use a method that begins with "distinct degree factorization", you only need the factor giving the product of the linear factors.
A: There is also the following way.
Let $k$ be an integer number.
Thus, $$x^4 + 21x^3 + 5x^2 + 7x + 1\equiv x^4-2x^3+5x^2+7x+1=$$
$$=(x^2-x+k)^2-((2k-4)x^2-(2k+7)x+k^2-1)).$$
Now, we'll choose a value of $k$ for which $$(2k+7)^2-8(k-2)(k^2-1)\equiv0.$$
We see that $k=6$ is valid.
Id est, $$x^4 + 21x^3 + 5x^2 + 7x + 1\equiv(x^2-x+6)^2-(8x^2-19x+35)\equiv$$
$$\equiv(x^2-x+6)^2-(100x^2-180x+81)=(x^2-x+6)^2-(10x-9)^2=$$
$$=(x^2-11x+15)(x^2+9x-3)$$ and the rest is smooth.
A: Let 
$$f(x)=x^4 -2x^3 + 5x^2 + 7x + 1\tag{1}$$
be defined over the finite field $\mathbb{F}_{23}$. Now check for a linear factor by checking for roots  over  $\mathbb{F}_{23}=\{0,\pm1,\pm2,\pm3,\pm4,\pm5\pm6,\pm7,\pm8,\pm9,\pm10,\pm11\}$. We find $f(-4)=f(-5)=0$, so $(x+4)$ and $(x+5)$ are linear factors. Now factor $f$ as two quadratics modulo $23$:
\begin{align*}
f(x)&=(x^2+9x-3)(x^2+ax+b)\\
&=x^4+(9+a)x^3+(9a-3+b)x^2+(9b-3a)x-3b
\end{align*}
Comparing the coefficients in $(1)$ for the powers of $x$:
\begin{array}\\
[x^3:] & -2=9+a\\
 [x^2:] & 5=9a-3+b\\
 [x:] & 7=9b-3a\\
 [const:]& 1=-3b\\
\end{array}
with $a$, $b$, $c$, $d\in\mathbb{F}_{23}$. Note this is a finite field so $-3b=1$ means $-3$ and $b$ are inverse mod $23$, making $b=15$. Now $a=-2-9=-11=12$ giving the factorization 
$$f(x)=(x^2+12x+15)(x+5)(x+4)$$
with the quadratic factor irreducible over $\mathbb{F}_{23}$ as it has no roots, since the discriminant of $(x^2+12x+15)$ is $15$ which is not a square modulo $23$.
