I did an implementation of the Tonelli-Shanks algorithm as defined on Wikipedia. I put it here for review and sharing purpose. [Legendre Symbol implementation][1]: def legendre_symbol(a, p): """ Legendre symbol Define if a is a quadratic residue modulo odd prime http://en.wikipedia.org/wiki/Legendre_symbol """ ls = pow(a, (p - 1)/2, p) if ls == p - 1: return -1 return ls [Prime modular square root][2] (I just renamed the solution variable R to x and n to a): def prime_mod_sqrt(a, p): """ Square root modulo prime number Solve the equation x^2 = a mod p and return list of x solution http://en.wikipedia.org/wiki/Tonelli-Shanks_algorithm """ a %= p # Simple case if a == 0: return [0] if p == 2: return [a] # Check solution existence on odd prime if legendre_symbol(a, p) != 1: return [] # Simple case if p % 4 == 3: x = pow(a, (p + 1)/4, p) return [x, p-x] # Factor p-1 on the form q * 2^s (with Q odd) q, s = p - 1, 0 while q % 2 == 0: s += 1 q //= 2 # Select a z which is a quadratic non resudue modulo p z = 1 while legendre_symbol(z, p) != -1: z += 1 c = pow(z, q, p) # Search for a solution x = pow(a, (q + 1)/2, p) t = pow(a, q, p) m = s while t != 1: # Find the lowest i such that t^(2^i) = 1 i, e = 0, 2 for i in xrange(1, m): if pow(t, e, p) == 1: break e *= 2 # Update next value to iterate b = pow(c, 2**(m - i - 1), p) x = (x * b) % p t = (t * b * b) % p c = (b * b) % p m = i return [x, p-x] If you have any optimization or found any error, please report it. [1]: http://en.wikipedia.org/wiki/Legendre_symbol [2]: http://en.wikipedia.org/wiki/Tonelli-Shanks_algorithm