I believe that the question whether $(p-1)! \bmod p^2$ is computable in polynomial time is still open. The best I am aware of is $O(p^{1/2+\epsilon})$, but I am not a specialist on this issue.

The computation of $(p-1)! \bmod p^2$ plays a role in certain p-adic algorithms to compute zeta functions of hypersurfaces in $\mathbb{P}^n$. There is an algorithm of David Harvey for zeta functions of hypersurfaces, where this question plays a arole. I.e., the fact that this is not known that $(p-1)! \bmod p^2$ can be computed in less than $O(p^{1/2+\epsilon})$ is an obstruction to improve the complexity of this algorithm. (See the slides "Computing zeta functions of projective surfaces in large characteristic" on http://www.cims.nyu.edu/~harvey/.)

I believe that the question whether $(p-1)! \bmod p^2$ is computable in polynomial time is still open. The best I am aware of is $O(p^{1/2+\epsilon})$, but I am not a specialist on this issue.

The computation of $(p-1)! \bmod p^2$ plays a role in certain p-adic algorithms to compute zeta functions of hypersurfaces in $\mathbb{P}^n$. There is an algorithm of David Harvey for zeta functions of hypersurfaces, where your question plays a role. I.e., the fact that this is not known that $(p-1)! \bmod p^2$ can be computed in less than $O(p^{1/2+\epsilon})$ is an obstruction to improve the complexity of this algorithm. (See the slides "Computing zeta functions of projective surfaces in large characteristic" on http://www.cims.nyu.edu/~harvey/.)