All square roots of $y$ are obtained by taking any particular square root of $y$ and multiply it by a square root of 1. So, the problem is split into two:
- find one (any) square root of $y$;
- find all square roots of 1.
In PARI/GP, both problems are solved more or less easily. E.g., for the first problem one can employ p-adic numbers and compute sqrt(y + O(p^t)).
This script computes all square roots, even modulo composite numbers.