math - Elliptic curve point addition over a finite field in Python -

in short, im trying add 2 points on elliptic curve y^2 = x^3 + ax + b on finite field fp. have working implementation on r, not know how alter general formulas ive found in order them sustain addition on fp.

when p not equal q, , z sum of p , q:

        dydx = (q.y - p.y)/(q.x - p.x)         z.x = dydx**2 - p.x - q.x         z.y = dydx * (z.x - p.x) + p.y 

when p equals q, again z sum:

        dydx = (3 * (p.x)**2 + self.a)/(2 * p.y)         z.x = dydx**2 - 2 * p.x         z.y = dydx * (z.x - p.x) + p.y 

these same formules found here. needs modified?

clarification: code above works elliptic curves on r. looking find needs changed work addition of points on finite field of order p.

there couple of issues here. first have wrong formulas: formulas negation of sum, or equivalently third point of curve lies on line through p , q. compare formula linked on wikipedia, , you'll see have z.y negation of value have.

a second issue formulas don't take origin account: if p inverse of q in group law, p + q origin, doesn't lie on affine part of curve , can't described pair of (x, y) coordinates. you'll need way specify origin.

let's write code. first need representation points on curve. use string 'origin' represent origin, , simple named tuple points on affine part of elliptic curve.

# create simple point class represent affine points. collections import namedtuple point = namedtuple("point", "x y")  # point @ infinity (origin group law). o = 'origin' 

for demonstration purposes, choose particular elliptic curve , prime. prime should greater 3 addition formulas valid. write function can use check particular point valid representation of point on curve. useful in checking didn't make mistakes in implementing addition formulas.

# choose particular curve , prime.  assume p > 3. p = 15733 = 1 b = 3  def valid(p):     """     determine whether have valid representation of point     on our curve.  assume x , y coordinates     reduced modulo p, can compare     2 points equality simple ==.     """     if p == o:         return true     else:         return (             (p.y**2 - (p.x**3 + a*p.x + b)) % p == 0 ,             0 <= p.x < p , 0 <= p.y < p) 

to divisions modulo p you'll need way compute inverses modulo p. simple , reasonably efficient trick here use python's three-argument variant of pow function: fermat's little theorem, pow(a, p-2, p) give inverse of a modulo p (provided of course that inverse exists - is, a not divisible p).

def inv_mod_p(x):     """     compute inverse x modulo p, assuming x     not divisible p.     """     if x % p == 0:         raise zerodivisionerror("impossible inverse")     return pow(x, p-2, p) 

and finally, here 2 functions compute negation , addition on elliptic curve. addition function based directly on formulas gave (after correcting sign of z.y), makes use of inv_mod_p perform divisions modulo p, , final reduction modulo p computed x , y coordinates.

def ec_inv(p):     """     inverse of point p on 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 points p , q on elliptic curve y^2 = x^3 + ax + b.     """     if not (valid(p) , valid(q)):         raise valueerror("invalid inputs")      # deal special cases either p, q, or p + q     # origin.     if p == o:         result = q     elif q == o:         result = p     elif q == ec_inv(p):         result = o     else:         # cases not involving 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)      # above computations *should* have given point     # on curve.     assert valid(result)     return result 

we can check code above creating handful of points on curve , checking obey 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 2 different ways. assert ec_add(p, threep) == ec_add(twop, twop) # check associative law. assert ec_add(p, ec_add(q, r)) == ec_add(ec_add(p, q), r)