  Wilson's theorem for the factorial got generalized to the moduli $p^2$ in
1900 and $p^3$ in 2000 by J.W.L. Glaisher and Z-H. Sun respectively. This paper
which studies more generally the multiple harmonic sums $\mathcal{H}_{\lbrace
s\rbrace^{2l}=1;p-1},2\leq 2l\leq p-1$ modulo $p^4$ in association with the
Stirling numbers $\left[\begin{array}{l}\;\;\;p\\2s-1\end{array}\right], 2\leq
2s\leq p-1$ modulo $p^4$ is concerned with establishing a generalization of
Wilson, Glaisher and Sun's results to the modulus $p^4$. We also break
p-residues of convolutions of three divided Bernoulli numbers of respective
orders $p-1$, $p-3$ and $p-5$ into smaller pieces and generalize some results
of Sun for some of the generalized harmonic numbers of order $p-1$ modulo
$p^4$.
