Implementation of the algorithm is in Sagemath:
x# Finding quadratic residue of a mod pfrom sage.rings.finite_rings.integer_mod import square_root_mod_prime a = <some a>p = <some p> print(square_root_mod_prime(Mod(a, p), p))
# Finding quadratic residue of a mod p
from sage.rings.finite_rings.integer_mod import square_root_mod_prime
a = <some a>
p = <some p>
print(square_root_mod_prime(Mod(a, p), p))