All square roots of $y$ are obtained by taking any particular square root of $y$ and multiplying 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](https://sagecell.sagemath.org/?z=eJx9lEuP0zAQgO-V-h8GVUh26642WThV4cSFQ1kJEJc8kJsmrLuxkzrukoL474zzdCqxPbTuvGc-j6MIeFFAfb5wnYEuS1NDmYMHsjyCn6jl4g8qtfHwv08UhQDklWi6Wy4ARE5UENwz0Jm5aEXCfXkk98yjMQXXwptbeMy_tfBHC2sQeuwhZu_oYNPLtuiYELX16GY8bb0Yj8ra_V0ulovotX6k2w2RfTOG5Uwz7bMvXTaD4hdeXLgRpSISi23F2kf5NAozlg8S60_M2MFLltalNqQQuSHap2MbOQbIeWpKLZRBtwi9epW2OTOrIpKbWvzOSE5D7E2wtvs2lk1OPNjAI8lDwbw4sT9-jBnYXAB92KGUITastL_2k5VmcGIgAnLCEb5F4Q4Owfcs1dkLOQjF9ZVYK0zVWkQr28YO2irSJ6GyOiM4kFBssEQYJX0WDP_MCPolh_A5Xmv8whHYMYyY_guqsaAuRQkSeA3mx6fP3_aPH2sgv4R5Qiz6KtTPzkZQi5PbqeA8LdFzT7RCFgyQSNGNofsgqjei7hKiKRuuZExdq2oG6XyHcV11YTP0hCpLyNGZCWJhwbW82plhYZZa1UGqempAb5x5VRVX0mw_WOQNa-6qpOJHkVY6S7v_SBramzc56iB0i8hLjTWQOnAqCXEnY8omI-uVlirleEFZOPBzPEwo4nWLsMYTxXV1Urjnr5npXoPXtk-NUEdiFpdiIOcPSl-bO4phK6LGomAnNq5DQzendSeVdi7uTdi34eXNjPEGrLR9sHr0ekZ-trfDdOh0Z6tubW_T4Kv0PnnoM7k2clKi7h_GXIXA&lang=gp&interacts=eJyLjgUAARUAuQ==) computes all square roots, even modulo composite numbers.