y^2 = x^3 + a*x + b
a = 2
b = 1
# pick a nice number from https://en.wikipedia.org/wiki/List_of_prime_numbers
p = 3571 # 941
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
y, x = np.ogrid[-5:5:100j, -5:5:100j]
plt.contour(x.ravel(), y.ravel(), pow(y, 2) - pow(x, 3) - x * a - b, [0])
plt.grid()
plt.show()
for i in range(0, p):
y2 = (pow(i, 3) + i * a + b) % p
y = sqrt(y2)
if y % 1 == 0:
y = prime_mod_sqrt(y2, p)
for yy in y:
#print '(%d, %d)' % (i, yy)
plt.plot(i, yy, 'ro')
#print '(Inf)'
#plt.axis([-1, 5, -1, 5])
plt.grid()
plt.show()
# via http://codereview.stackexchange.com/questions/43210/tonelli-shanks-algorithm-implementation-of-prime-modular-square-root
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
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]