Elliptical Curve point addition and doubling on montgomery curve

Question

I am trying to teach myself about elliptical curve cryptography using online resources, and have created a program to perform point addition and point doubling on a montgomery curve. To test my function, I am calculating 4p using two methods: 2(2p) and p + p + p + p. To my understanding, the values should be equivalent, but I'm getting different values.

I'm using curve M-221 defined as:

y^2 = x^3 + 117050 * x^2 + x with field 2^221 - 3 and the generator point:

{x: 1868095183606742538701435617698957211251810735713795112730238517387, y: 1432401617868315471469241676678456314840808294315547381535466705153}

When running my test I get:

p4_double:   {'x': 1158296677984570086957317285562479367861170755348925863907217477275, 'y': 1551028763014061658160471452226999211325397814089364066655644963406}
p4_add:      {'x': 27500914413453982607377956096660382032281000528958526135014696542, 'y': 1168146124709862201096710659924046078056321465610616206479606975751}
--------------------------------------------------
p2:          {'x': 339504813444055524776350731343982961171002195927973873104309035752, 'y': 2932526515499324705932531354096697705273063568784183338538040616672}
p3:          {'x': 394721346804099203910004476506708200743020785679182057269678209841, 'y': 1663472166528890517288736425441643076121539170204221488121613770082}

My question is, what am I doing to make 2(2p) != p + p + p + p. I believe that my problem lies in one of the following:

• Incorrect mathematics (which is why I believe this question belongs here and not in crypto)
• use of rounding
• use of modulo

Interestingly, 2p and 4p_add are very similar, lending me to believe the error lies in adding two points. I appreciate all feedback, as I am very new to elliptical curves, and especially with working within finite fields.

Edit: So I found out that I needed to use modulo much more frequently throughout my program. I now have a function that works on small fields, but still fails the given curve. Given these parameters the function works perfectly - even having 16p_add and 16p_double coinciding:

curve = {
    "A": 5,
    "B": 3,
    "field": 65537,
    "generator": {
        "x": 3,
        "y": 5
    }
}

I therefore believe that the error is coming from python not keeping enough significant bits on large numbers. However, I would still love help to fix my program if the error lies elsewhere. (At this point this question may belong in a section dedicated to python/programming) Here is my code:

# Provides elliptical-curve math capability for deffie-hellman key exchange

# The initial parameters for the montgomery curve: B * y^2 = x^3 + A * X^2 + x
# values obtained for curve M-221 (2013 Aranha–Barreto–Pereira–Ricardini)
# Generator point obtained from tutorial at: https://sage.math.clemson.edu:34567/home/pub/847/
curve = {
    "A": 117050,
    "B": 1,
    "field": 2 ** 221 - 3,
    "generator": {
        "x": 1868095183606742538701435617698957211251810735713795112730238517387,
        "y": 1432401617868315471469241676678456314840808294315547381535466705153
    }

}


def gcd_extended(a, b):
    """
    Uses Euclidean algorithm
    :param a: Int the value to test
    :param b: Int the modulo to use
    :return: x, y
    """
    if a == 0:
        return 0, 1
    x, y = gcd_extended(b % a, a)

    return y - (b // a) * x, x


def inverse_mod(a, b):
    """
    Calculates the inverse mod for a given a^-1 % b
    :param a: Int the value to find the inverse mod of
    :param b: Int the modulo to use
    :return: Int the inverse modulo
    """
    x, y = gcd_extended(a, b)
    # final modulo is to ensure x lies within the modulo field
    return x % b


def point_combination(_m, _p1, _p2):
    """
    Calculates p1 + p2 = p3 returning p3 - used by both doubling and adding functions
    I separated this function from the other two in hopes to make the code more readable
    :param _m: The slope of the line through p1, p2 or tangent line of p1
    :param _p1: Dictionary the first point given as {"x", "y"}
    :param _p2: Dictionary the first point given as {"x", "y"} or p1 if doubling
    :return: p3 as Dictionary {"x", "y"} within the curve field
    """

    # Calculate x, note that _p2["x"] == _p1["x"] if we are doubling
    # I round here to keep x as an int and take modulo field at this point
    x = (curve["B"] * _m ** 2 - (_p1["x"] + _p2["x"]) - curve["A"]) \
        % curve["field"]

    # Calculate y from x
    # I round here to keep y as an int and take modulo field at this point
    y = (_m * (_p1["x"] - x) - _p1["y"]) % curve["field"]

    # return p3 - note that initially we have R, so flip y along x axis for p3
    return {"x": x, "y": y}


def points_add(_p1, _p2):
    """
    Calculates p1 + p2 = p3 returning p3
    :param _p1: Dictionary the first point given as {"x", "y"}
    :param _p2: Dictionary the second point given as {"x", "y"}
    :return: p3 as Dictionary {"x", "y"} within the curve field
    """

    # Calculate the slope through points _p1 and _p2
    _m = (_p2["y"] - _p1["y"]) * inverse_mod(_p2["x"] - _p1["x"], curve["field"])

    # use slope and two points to find p3 and return it
    return point_combination(_m, _p1, _p2)


def point_double(_p):
    """
    Calculates p + p = 2p returning 2p
    :param _p: Dictionary the point to double given as {"x", "y"}
    :return: 2p as Dictionary {"x", "y"} within the curve field
    """

    # Calculate the tangent line using the derivative dx/dy
    _m = ((3 * _p["x"] ** 2 + 2 * curve["A"] * _p["x"] + 1) *
          inverse_mod(2 * curve["B"] * _p["y"], curve["field"])) % curve["field"]

    # use slope and the point to find p3 and return it
    # p is listed twice in this function, as the formula for x uses the x value of two points
    # in point doubling the second x value = p["x"]
    return point_combination(_m, _p, _p)


# just shorthand for the generator point
p = curve["generator"]

# calculate 4p via point doubling twice and adding p 4 times


# start with an initial doubling
p2 = point_double(p)

p4_double = point_double(p2)

p3 = points_add(p, p2)
p4_add = points_add(p, p3)

# in theory the values should be the same as 2 * 2p  ==  p + p + p + p = 2p + p + p
print("p4_double: \t", p4_double)
print("p4_add: \t", p4_add)

print("-" * 50)
print("p2: \t\t", p2)
print("p3: \t\t", p3)

Answer 1

As I mentioned in my edit, I was able to successfully edit my function to work without any problems on a smaller curve, using the parameters:

curve = {
    "A": 5,
    "B": 3,
    "field": 65537,
    "generator": {
        "x": 3,
        "y": 5
    }
}

The two functions add and double are consistent to 2^14p, that is 14 doublings. After this point, the two values begin to deviate due to significant bit errors. I tried changing B -> 1 and found the same error occurs after 2^14p.

The initial problem was that I needed to find the inverse modulo whenever I divided, instead of rounding after division.

At this point, there appears to be nothing wrong with the math. Here is the final code for those who are curious:

# Provides elliptical-curve math capability for deffie-hellman key exchange

# The initial parameters for the montgomery curve: B * y^2 = x^3 + A * X^2 + x
# values obtained for curve M-221 (2013 Aranha–Barreto–Pereira–Ricardini)
# Generator point obtained from tutorial at: https://sage.math.clemson.edu:34567/home/pub/847/
# These particular values where taken from: 
# https://math.stackexchange.com/questions/1396544/problem-with-elliptic-curve-in-montgomery-form?rq=1
curve = {
    "A": 5,
    "B": 1,
    "field": 65537,
    "generator": {
        "x": 3,
        "y": 5
    }
}


def gcd_extended(a, b):
    """
    Uses Euclidean algorithm
    :param a: Int the value to test
    :param b: Int the modulo to use
    :return: x, y
    """
    if a == 0:
        return 0, 1
    x, y = gcd_extended(b % a, a)

    return y - (b // a) * x, x


def inverse_mod(a, b):
    """
    Calculates the inverse mod for a given a^-1 % b
    :param a: Int the value to find the inverse mod of
    :param b: Int the modulo to use
    :return: Int the inverse modulo
    """
    x, y = gcd_extended(a, b)
    # final modulo is to ensure x lies within the modulo field
    return x % b


def point_combination(_m, _p1, _p2):
    """
    Calculates p1 + p2 = p3 returning p3 - used by both doubling and adding functions
    I separated this function from the other two in hopes to make the code more readable
    :param _m: The slope of the line through p1, p2 or tangent line of p1
    :param _p1: Dictionary the first point given as {"x", "y"}
    :param _p2: Dictionary the first point given as {"x", "y"} or p1 if doubling
    :return: p3 as Dictionary {"x", "y"} within the curve field
    """

    # Calculate x, note that _p2["x"] == _p1["x"] if we are doubling
    x = (curve["B"] * _m ** 2 - (_p1["x"] + _p2["x"]) - curve["A"]) \
        % curve["field"]

    # Calculate y from x
    y = (_m * (_p1["x"] - x) - _p1["y"]) % curve["field"]

    return {"x": x, "y": y}


def points_add(_p1, _p2):
    """
    Calculates p1 + p2 = p3 returning p3
    :param _p1: Dictionary the first point given as {"x", "y"}
    :param _p2: Dictionary the second point given as {"x", "y"}
    :return: p3 as Dictionary {"x", "y"} within the curve field
    """

    # Calculate the slope through points _p1 and _p2
    # This is one location where inverese_mod was needed as we cannot simply find
    # (y2-y1)/(x2-x1) - since this is % field we need to multiply by inverse mod (x2-x1) 
    _m = (_p2["y"] - _p1["y"]) * inverse_mod(_p2["x"] - _p1["x"], curve["field"])

    # use slope and two points to find p3 and return it
    return point_combination(_m, _p1, _p2)


def point_double(_p):
    """
    Calculates p + p = 2p returning 2p
    :param _p: Dictionary the point to double given as {"x", "y"}
    :return: 2p as Dictionary {"x", "y"} within the curve field
    """

    # Calculate the tangent line using the derivative dy/dx
    # Another location where we need to use inverse_mod as dy/dx % field ->
    # dy * inverse_mod(dx)
    _m = ((3 * _p["x"] ** 2 + 2 * curve["A"] * _p["x"] + 1) *
          inverse_mod(2 * curve["B"] * _p["y"], curve["field"])) % curve["field"]

    # use slope and the point to find p3 and return it
    # p is listed twice in this function, as the formula for x uses the x value of two points
    # in point doubling the second x value = p["x"]
    return point_combination(_m, _p, _p)


# just shorthand for the generator point
p = curve["generator"]

# calculate 4p via point doubling twice and adding p 4 times


# start with an initial doubling
p2 = point_double(p)
# 2^n will be the value (2^n)P to find
n = 14

# For the doubling function
pn_pow = p2
for i in range(n-1):
    pn_pow = point_double(pn_pow)
    print(f"p{2 ** (i + 1)}: \t\t", pn_pow)

# for the adding
pn_add = p2
for i in range(2 ** n - 2):
    pn_add = points_add(p, pn_add)

# the values are equivalent to (2^14)p and deviate by (2^15)p
print(f"p{2 ** n}_double: \t", pn_pow)
print(f"p{2 ** n}_add:    \t", pn_add)

``` 

Answer 2

It may be a good idea to use already written code. Sage is the canonical choice, elliptic curve operations are supported in an object oriented way, the code is python "with batteries" (finally performing operations in other computer algebra systems like pari/gp).

Here is all what is needed to compute $4P$ either as $P+P+P+P$, or as $2(2P)$ in the given example.

To initialize a curve defined by $y^2 +a_1xy +a_3 y = x^3 + a_2x^2+a_4x+a_6$ just pass the list $[a_1,a_2,a_3,a_4,a_6]$ to the constructor as below. (If the curve has $a_1=a_3=a_2=0$, then one can pass simply $[s_4,a_6]$ instead. This is not the case for our curve $M-221$, but it is good to know.) The code and the computation, as shown in the IPython console is as follows...

---

F = GF(2^221 - 3) 
E = EllipticCurve(F, [0, 117050, 0, 1, 0])
x0 = 1868095183606742538701435617698957211251810735713795112730238517387
y0 = 1432401617868315471469241676678456314840808294315547381535466705153
P = E.point( (x0, y0) )    # and the point is also checked
P + P + P + P == 2*(2*P)

The last line delivers...

sage: P + P + P + P == 2*2*P
True
sage: P + P + P + P == 4*P
True

And the point is

sage: Q = 4*P    # its components are - in case Q is not (0:1:0)
sage: Q.xy()[0]  # first component, x component of Q
1158296677984570086957317285562479367861170755348925863907217477275
sage: Q.xy()[1]  # second component, y component of Q
1551028763014061658160471452226999211325397814089364066655644963406

Since the library is free, and the whole involved code is accessible, it may be a good idea to just import it...