I'm trying to implement a simple elliptic curve encryption program but I can't get the expected output of doubling and adding a Point P till 12P .The curve equation isy^2 = x^3 +ax + b mod p. According to this site 3P = [10, 6] when P = [5, 1] while I get 3p = [10, 5]. The equations I use can be found on Wikipedia.
P = [5, 1]
prime = 17
a = 2
b = 2
def gcdExtended(a, b):
if a == 0:
return b, 0, 1
gcd, x1, y1 = gcdExtended(b % a, a)
x = y1 - (b // a) * x1
y = x1
return gcd, x, y
def double_point(point: list):
x = point[0]
y = point[1]
s = ((3*(x**2)+a) * (gcdExtended(2*y, prime)[1])) % prime
newx = (s**2 - x - x) % prime
newy = (s * (x - newx) - y) % prime
return [newx, newy]
def add_points(P: list, Q: list):
x1 = P[0]
y1 = P[1]
x2 = Q[0]
y2 = Q[1]
s = ((y2 - y1) * ((gcdExtended(x2-x1, prime))[1] % prime)) % prime
newx = (s**2 - x1 - x2) % prime
newy = (s * (x1 - newx) - y1) % prime
return [newx, newy]
Q = P
index = 2
while True:
if Q[0] == P[0] and Q[1] == P[1]:
print("doubling")
Q = double_point(P)
else:
print("adding")
Q = add_points(Q, P)
if index == 12 :
break
print(f"{index}P = {Q}")
index += 1
If the point [5,1] is added successively, the following sequence is obtained:
1P = [ 5, 1]
2P = [ 6, 3]
3P = [10, 6]
4P = [ 3, 1]
5P = [ 9, 16]
6P = [16, 13]
7P = [ 0, 6]
8P = [13, 7]
9P = [ 7, 6]
10P = [ 7, 11]
11P = [13, 10]
12P = [ 0, 11]
13P = [16, 4]
14P = [ 9, 1]
15P = [ 3, 16]
16P = [10, 11]
17P = [ 6, 14]
18P = [ 5, 16]
19P = point at infinity
This can be verified e.g. here.
The problem in the posted code is that the method to determine the modular inverse, gcdExtended(a, b), is only valid for positive a and b. While in double_point and add_points b has the value prime (= 17 > 0), a can take negative values.
gcdExtended generally returns wrong values for negative a:
gcdExtended returns for these values: gcdExtended(5, 17)[1] = 7 (which is true) and gcdExtended(-12, 17)[1] = -7 (which is false).To allow negative values for a, e.g. the following methods can be defined, see here:
def sign(x):
return 1 if x >= 0 else -1
def gcdExtendedGeneralized(a, b):
gcd, x1, y1 = gcdExtended(abs(a), b)
return gcd, (sign(a) * x1) % b, y1 % b
Replacing gcdExtended with gcdExtendedGeneralized in double_point and add_points provides the correct values (note that the current implementation does not consider the point at infinity).