There are a couple of issues here. First is that you have the wrong formulas: those are the formulas for the negation of the sum, or equivalently the third point of the curve that lies on the line through P and Q. Compare with the formula you linked to on Wikipedia, and you'll see that what you have for Z.y is the negation of the value they have.
A second issue is that your formulas don't take the origin into account: if P is the inverse of Q in the group law, then P + Q will be the origin, which doesn't lie on the affine part of the curve and so can't be described as a pair of (x, y) coordinates. So you'll need some way to specify the origin.
Let's write some code. First we need a representation for the points on the curve. We use the string 'Origin' to represent the origin, and a simple named tuple for the points on the affine part of the elliptic curve.
# Create a simple Point class to represent the affine points.
from collections import namedtuple
Point = namedtuple("Point", "x y")
# The point at infinity (origin for the group law).
O = 'Origin'
For demonstration purposes, we choose a particular elliptic curve and prime. The prime should be greater than 3 for the addition formulas to be valid. We also write a function that we can use to check that a particular point is a valid representation of a point on the curve. This will be useful in checking that we didn't make any mistakes in implementing the addition formulas.
# Choose a particular curve and prime. We assume that p > 3.
p = 15733
a = 1
b = 3
def valid(P):
"""
Determine whether we have a valid representation of a point
on our curve. We assume that the x and y coordinates
are always reduced modulo p, so that we can compare
two points for equality with a simple ==.
"""
if P == O:
return True
else:
return (
(P.y**2 - (P.x**3 + a*P.x + b)) % p == 0 and
0 <= P.x < p and 0 <= P.y < p)
To do divisions modulo p you'll need some way to compute inverses modulo p. A simple and reasonably efficient trick here is to use Python's three-argument variant of the pow function: by Fermat's Little Theorem, pow(a, p-2, p) will give an inverse of a modulo p (provided of course that that inverse exists - that is, that a is not divisible by p).
def inv_mod_p(x):
"""
Compute an inverse for x modulo p, assuming that x
is not divisible by p.
"""
if x % p == 0:
raise ZeroDivisionError("Impossible inverse")
return pow(x, p-2, p)
And finally, here are the two functions to compute negation and addition on the elliptic curve. The addition function is based directly on the formulas you gave (after correcting the sign of Z.y), makes use of inv_mod_p to perform the divisions modulo p, and does a final reduction modulo p for the computed x and y coordinates.
def ec_inv(P):
"""
Inverse of the point P on the elliptic curve y^2 = x^3 + ax + b.
"""
if P == O:
return P
return Point(P.x, (-P.y)%p)
def ec_add(P, Q):
"""
Sum of the points P and Q on the elliptic curve y^2 = x^3 + ax + b.
"""
if not (valid(P) and valid(Q)):
raise ValueError("Invalid inputs")
# Deal with the special cases where either P, Q, or P + Q is
# the origin.
if P == O:
result = Q
elif Q == O:
result = P
elif Q == ec_inv(P):
result = O
else:
# Cases not involving the origin.
if P == Q:
dydx = (3 * P.x**2 + a) * inv_mod_p(2 * P.y)
else:
dydx = (Q.y - P.y) * inv_mod_p(Q.x - P.x)
x = (dydx**2 - P.x - Q.x) % p
y = (dydx * (P.x - x) - P.y) % p
result = Point(x, y)
# The above computations *should* have given us another point
# on the curve.
assert valid(result)
return result
We can check the code above by creating a handful of points on the curve and checking that they obey the expected arithmetic laws.
P = Point(6, 15)
Q = Point(8, 1267)
R = Point(2, 3103)
TwoP = ec_add(P, P)
ThreeP = ec_add(TwoP, P)
# Compute 4P two different ways.
assert ec_add(P, ThreeP) == ec_add(TwoP, TwoP)
# Check the associative law.
assert ec_add(P, ec_add(Q, R)) == ec_add(ec_add(P, Q), R)
