Elliptic curve double and add implementation in python

Question

I am trying to implement the "double and add" algorithm to quickly multiply points on an elliptic curve in Python (3, please).

Based off this previous answer (about addition and doubling),

Elliptic curve point addition over a finite field in Python

the Wikipedia page,

https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication

and my textbook (Information Security, by Mark Stamp), I came up with the following code:

def point_add(N_x, N_y, Q_x, Q_y, p):
    m = (Q_y - N_y) * pow((Q_x-N_x), p-2, p)
    ret_x = (m ** 2 - N_x - Q_x) % p
    ret_y = (m*(N_x - ret_x) - N_y) % p
    return ret_x, ret_y


def point_double(N_x, N_y, a, p):
    m = (3*(N_x ** 2)+a) * pow(2*N_y, p-2, p)   
    ret_x = (m ** 2 - N_x - Q_x) % p
    ret_y = (m*(N_x - ret_x) - N_y) % p
    return ret_x, ret_y


multiplier_A = 44
multiplier_B = 57
a = 10
b = -21
p = 41
N_x = 3
N_y = 6
Q_x = 0
Q_y = 0

multiplier_A = list(bin(multiplier_A)[2:])
multiplier_B = list(bin(multiplier_B)[2:])

for x_a in multiplier_A:
    if x_a == '1':
        Q_x, Q_y = point_add(N_x, N_y, Q_x, Q_y, p)
    N_x, N_y = point_double(N_x, N_y, a, p)

print(Q_x, Q_y)

Which sure enough returns a wrong result.

Where do I go wrong? A new implementation from scratch is also more than welcome.

UPDATE

I want to multiply the point (3,6) times multiplier_A and multiplier_B. For now I'm trying to only do multiplier_A(3,6).

a, b, and p are the elliptic curve paramethers

Answer 1

Both of your functions are totally correct, except small typo that point_double() should have N_x instead of Q_x inside its body, because there is no Q_x argument in function (Q_x is undefined).

I checked them against my own implementation that I used many times in different algorithms, my implementation is called add_dbl() in code below.

You can see below that I used standard Secp256k1 elliptic curve, that is default curve for BitCoin.

Function test() generates 200 random points and does addition and doubling with both my and yours implementations, and checks if they give equal results.

Code doesn't have any non-standard dependecies, just imports standard random module.

Try it online!

def get_y(A, B, N, x):
    assert N % 4 == 3, N
    y2 = x ** 3 + A * x + B
    y = pow(y2, (N + 1) // 4, N)
    if pow(y, 2, N) != y2 % N:
        return None
    return y

def add_dbl(A, N, Px, Py, Qx = None, Qy = None):
    if Qx is None and Qy is None:
        Qx, Qy = Px, Py
    if Px == Qx and Py == Qy:
        s = ((Px * Px * 3 + A) * pow(Py * 2, -1, N)) % N
    else:
        s = ((Py - Qy) * pow(Px - Qx, -1, N)) % N
    x = (s * s - Px - Qx) % N
    y = (s * (Px - x) - Py) % N
    return x, y

def base_point():
    # https://en.bitcoin.it/wiki/Secp256k1
    # https://www.secg.org/sec2-v2.pdf
    A = 0
    B = 7
    N = 0xFFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFE_FFFFFC2F
    x = 0x79BE667E_F9DCBBAC_55A06295_CE870B07_029BFCDB_2DCE28D9_59F2815B_16F81798
    y = 0x483ADA77_26A3C465_5DA4FBFC_0E1108A8_FD17B448_A6855419_9C47D08F_FB10D4B8
    q = 0xFFFFFFFF_FFFFFFFF_FFFFFFFF_FFFFFFFE_BAAEDCE6_AF48A03B_BFD25E8C_D0364141
    return A, B, N, x, y

def point_add(N_x, N_y, Q_x, Q_y, p):
    m = (Q_y - N_y) * pow((Q_x-N_x), p-2, p)
    ret_x = (m ** 2 - N_x - Q_x) % p
    ret_y = (m*(N_x - ret_x) - N_y) % p
    return ret_x, ret_y

def point_double(N_x, N_y, a, p):
    m = (3*(N_x ** 2)+a) * pow(2*N_y, p-2, p)   
    ret_x = (m ** 2 - N_x - N_x) % p
    ret_y = (m*(N_x - ret_x) - N_y) % p
    return ret_x, ret_y

def test():
    import random
    pb = base_point()
    A, B, N, *_ = pb
    for itest in range(200):
        while True:
            x0, x1 = [random.randrange(1, N) for i in range(2)]
            y0, y1 = [get_y(A, B, N, x) for x in (x0, x1)]
            if y0 is not None and y1 is not None:
                break
        
        add_a = add_dbl(A, N, x0, y0, x1, y1)
        add_b = point_add(x0, y0, x1, y1, N)
        assert add_a == add_b
        
        dbl_a = add_dbl(A, N, x0, y0)
        dbl_b = point_double(x0, y0, A, N)
        assert dbl_a == dbl_b

if __name__ == '__main__':
    test()