97 lines
2.3 KiB
Python
97 lines
2.3 KiB
Python
|
# from http://eli.thegreenplace.net/2009/03/07/computing-modular-square-roots-in-python/
|
||
|
|
|
||
|
|
|
||
|
|
def modular_sqrt(a, p):
|
||
|
|
""" Find a quadratic residue (mod p) of 'a'. p
|
||
|
|
must be an odd prime.
|
||
|
|
|
||
|
|
Solve the congruence of the form:
|
||
|
|
x^2 = a (mod p)
|
||
|
|
And returns x. Note that p - x is also a root.
|
||
|
|
|
||
|
|
0 is returned is no square root exists for
|
||
|
|
these a and p.
|
||
|
|
|
||
|
|
The Tonelli-Shanks algorithm is used (except
|
||
|
|
for some simple cases in which the solution
|
||
|
|
is known from an identity). This algorithm
|
||
|
|
runs in polynomial time (unless the
|
||
|
|
generalized Riemann hypothesis is false).
|
||
|
|
"""
|
||
|
|
# Simple cases
|
||
|
|
#
|
||
|
|
if legendre_symbol(a, p) != 1:
|
||
|
|
return 0
|
||
|
|
elif a == 0:
|
||
|
|
return 0
|
||
|
|
elif p == 2:
|
||
|
|
return p
|
||
|
|
elif p % 4 == 3:
|
||
|
|
return pow(a, (p + 1) / 4, p)
|
||
|
|
|
||
|
|
# Partition p-1 to s * 2^e for an odd s (i.e.
|
||
|
|
# reduce all the powers of 2 from p-1)
|
||
|
|
#
|
||
|
|
s = p - 1
|
||
|
|
e = 0
|
||
|
|
while s % 2 == 0:
|
||
|
|
s /= 2
|
||
|
|
e += 1
|
||
|
|
|
||
|
|
# Find some 'n' with a legendre symbol n|p = -1.
|
||
|
|
# Shouldn't take long.
|
||
|
|
#
|
||
|
|
n = 2
|
||
|
|
while legendre_symbol(n, p) != -1:
|
||
|
|
n += 1
|
||
|
|
|
||
|
|
# Here be dragons!
|
||
|
|
# Read the paper "Square roots from 1; 24, 51,
|
||
|
|
# 10 to Dan Shanks" by Ezra Brown for more
|
||
|
|
# information
|
||
|
|
#
|
||
|
|
|
||
|
|
# x is a guess of the square root that gets better
|
||
|
|
# with each iteration.
|
||
|
|
# b is the "fudge factor" - by how much we're off
|
||
|
|
# with the guess. The invariant x^2 = ab (mod p)
|
||
|
|
# is maintained throughout the loop.
|
||
|
|
# g is used for successive powers of n to update
|
||
|
|
# both a and b
|
||
|
|
# r is the exponent - decreases with each update
|
||
|
|
#
|
||
|
|
x = pow(a, (s + 1) / 2, p)
|
||
|
|
b = pow(a, s, p)
|
||
|
|
g = pow(n, s, p)
|
||
|
|
r = e
|
||
|
|
|
||
|
|
while True:
|
||
|
|
t = b
|
||
|
|
m = 0
|
||
|
|
for m in xrange(r):
|
||
|
|
if t == 1:
|
||
|
|
break
|
||
|
|
t = pow(t, 2, p)
|
||
|
|
|
||
|
|
if m == 0:
|
||
|
|
return x
|
||
|
|
|
||
|
|
gs = pow(g, 2 ** (r - m - 1), p)
|
||
|
|
g = (gs * gs) % p
|
||
|
|
x = (x * gs) % p
|
||
|
|
b = (b * g) % p
|
||
|
|
r = m
|
||
|
|
|
||
|
|
|
||
|
|
def legendre_symbol(a, p):
|
||
|
|
""" Compute the Legendre symbol a|p using
|
||
|
|
Euler's criterion. p is a prime, a is
|
||
|
|
relatively prime to p (if p divides
|
||
|
|
a, then a|p = 0)
|
||
|
|
|
||
|
|
Returns 1 if a has a square root modulo
|
||
|
|
p, -1 otherwise.
|
||
|
|
"""
|
||
|
|
ls = pow(a, (p - 1) / 2, p)
|
||
|
|
return -1 if ls == p - 1 else ls
|