As you might already know, solutions to the quintic can be expressed in terms of either ${}_4 F_3$ hypergeometric functions or Jacobi theta functions. See King or Prasolov/Solovyev for details.
For polynomials of higher degree, there is also a general formula for the roots, due to Umemura. The formulae involve the multidimensional generalization of the Jacobi theta functions (the Riemann theta function), and are a bit unwieldy; see Umemura's paper if you want more details. See also this preprint for a solution of the reduced polynomial equation $x^n-x-\alpha=0$ in terms of hypergeometric functions.
Let's start off by looking at the characteristic polynomial of the $2 \times 2$ matrix
$A = \begin{bmatrix} 0 & 1 \\ -a & -b \end{bmatrix}: \tag 1$
$\det (A - \lambda I) = \det \left (\begin{bmatrix} -\lambda & 1 \\ -a & -b - \lambda \end{bmatrix} \right ) = \lambda^2 + b\lambda + a; \tag 2$
we see from (2) that we may always present a $2 \times 2$ matrix with given characteristic polynomial $\lambda^2 + b\lambda + a$ in the form $A$; for example, if he quadratic is $\lambda^2 + 5\lambda + 1$, as in the present problem (tho' I have replaced $p$ with $\lambda$), we may take
$P = \begin{bmatrix} 0 & 1 \\ -1 & -5 \end{bmatrix}, \tag 3$
which may be easily checked:
$\det(P - \lambda I) = -\lambda(-5 - \lambda) - 1 ( -1) = \lambda^2 + 5\lambda + 1. \tag 4$
The matrices $A$ and $P$ above are part of a general paradigm for constructing matrices with a given characteristic polynomial, and it extends to higher dimensions. Now, there are a very many matrices possessed of a given characteristic polynomial, since it is a similarity invariant; that is, the characteristic polynomials of $X$ and $S^{-1}XS$ are always the same; thus it behooves us to find a matrix of particularly simple, general form for a given polynomial. If
$q(\lambda) = \displaystyle \sum_1^n q_i \lambda^i, \; q_n = 1, \tag5$
we define $C(q(\lambda))$ to be the $n \times n$ matrix
$C(q(\lambda)) = \begin{bmatrix} 0 & 1 & 0 & \ldots & 0 \\ 0 & 0 & 1 & \ldots & 0 \\
\vdots & \vdots & \ldots & \vdots & \vdots \\ -q_0 & -q_1 & -q_2 & \ldots & -q_{n - 1} \end{bmatrix}; \tag 6$
that is, $C(q(\lambda))$ has all $1$s on the superdiagonal, the negatives of the coefficients of $q(\lambda)$ on the $n$-th row, and $0$s everywhere else. We have
$C(q(\lambda) - \lambda I = \begin{bmatrix} -\lambda & 1 & 0 & \ldots & 0 \\ 0 & -\lambda & 1 & \ldots & 0 \\
\vdots & \vdots & \ldots & \vdots & \vdots \\ -q_0 & -q_1 & -q_2 & \ldots & -q_{n - 1} - \lambda \end{bmatrix}; \tag 7$
it is easy to see, by expanding in minors along the $n$-th row, that
$\det(C(q(\lambda)) - \lambda I) = q(\lambda); \tag 8$
also, since a matrix and its transpose have equal determinants, the transposed form $C(q(\lambda))$, $C^T(q(\lambda))$, also gives rise to the same characteristic polynomial. The matrices $C(q(\lambda))$ and $C^T(q(\lambda))$ are known as companion matrices for the polynomial $q(\lambda)$; the linked article has more of the story.
Best Answer
If a matrix $A\in M_n(\mathbb{Q})$ is random, then roughky speaking, its characteristic polynomial is random. Then we may consider a (non monic) random polynomial $p=\sum_{0\leq i\leq 5}a_ix^i\in \mathbb{Z}_5[x]$.
The simplest method is to consider a positive integer $n$ and to randomly choose (independently) the $(a_i)$ uniformly in $\{-n..n\}$. Let $P_n$ be the associated probability that $p$ is irreducible and has $S_5$ as Galois group.
EDIT. Then $\lim_{n\rightarrow +\infty}P_n=1$. About this result, you can read
[1] J.P. Serre: Topics in Galois Theory.(the reading is hard)
[2] Igor Irvine: Galois groups of generic polynomials.
https://arxiv.org/pdf/1511.06446.pdf
A difficult problem is to estimate the speed of convergence towards $1$ of $P_n$; an upper bound is given in [1] and more precisely in [2].
To give an idea, here are the results of some random tests
$P_1\approx 28$%,$P_{10}\approx 85$%$,P_{100}\approx 98.2$%,$P_{1000}\approx 99.79$%.
Of course, $P_n$ depends on the degree of the polynomial $p$. When the degree increases, $P_n$ increases too.