How to add elliptic curve points in python?

Question

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

Answer 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:

• The modular inverse of 5 or -12 is 7: 5 x 7 mod17 = 35 mod17 = 1 and 7 x (-12) mod17 = -84 mod17 = 85 mod17 = 1.
• The 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).

Answer 2

You interchanged P and Q in add_points. Also a small simplification in your calculation of s:

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
    s = (y2-y1) * (gcdExtended(x2-x1, prime)[1] % 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(P, Q)

    if index == 12 :
        break

    print(f"{index}P = {Q}")
    index += 1

which results in

doubling
2P = [6, 3]
adding
3P = [10, 6]
adding
4P = [3, 1]
adding
5P = [9, 16]
adding
6P = [16, 13]
adding
7P = [0, 6]
adding
8P = [13, 8]
adding
9P = [8, 7]
adding
10P = [8, 10]
adding
11P = [13, 9]
adding