Multiple harmonic sums $\mathcal{H}_{\lbrace s\rbrace^{2l}=1;p-1}$ modulo $p^4$ and applications

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$.