Wilson's theorem modulo p² derived from Faulhaber polynomials

First, we present a new proof of Glaisher's formula dating from 1900 and concerning Wilson's theorem modulo p^2. Our proof uses p-adic numbers and Faulhaber's formula for the sums of powers (17th century), as well as more recent results on Faulhaber's coefficients obtained by Gessel and Viennot. Second, by using our method, we find a simpler proof than Sun's proof regarding a formula for (p-1)! modulo p^3, and one that can be generalized to higher powers of p. Third, we can derive from our method a way to compute the Stirling numbers modulo p^3, thus improving Glaisher and Sun's own results from 120 years ago and 20 years ago respectively. Last, our method allows to find new congruences on convolution of divided Bernoulli numbers and convolutions of divided Bernoulli numbers with Bernoulli numbers.