Implementing Elliptic Curve Cryptography

CS 463/480 Lecture, Dr. Lawlor

Our bottom line from the math of elliptic curves is to add points (x1,y1) to (x2,y2) on the elliptic curve y^2 = x^3 + a x + b, we compute a slope m and offset v:

Different Points
m = dy/dx = (y1-y2) / (x1-x2)

Same Point (tangent)
m = dy/dx = (3 x ^2 + a) / (2 y)

v = y1 - m*x1
x3 = m^2 - x1 - x2
y3 = - (m x3 + v)

The big complication involves the division:

All the arithmetic is modulo a big prime number p, including the division, which we need to implement via a modular inverse, computed via the extended euclidean algorithm.
If you would be dividing by zero, there are three possibilities:

(x1-x2)==0 and (y1-y2)==0, because you're adding a point to itself. Switch to the tangent code and continue.
(x1-x2)==0, but (y1-y2)!=0. These points are additive inverses, so their sum is the special "identity" point. I usually represent the identity as x=p, y=0.
(2y)==0. For most elliptic curves, this can only happen if you're adding the identity to itself.

One annoyance is in the standard C implementation of the % operator, an expression like (-(m x3 + v))%p gives a negative result, but modular arithmetic is never negative. This means we often need to add extra copies of p so the numbers always stay positive. Python is among the few languages where % always returns a positive result, but I've left the extra copies of p for portability to other languages.

Here's some python3 code to directly implement elliptic curve point addition and multiplication, including the special cases with the identity element:

# An elliptic curve has these fields:
#   p: the prime used to mod all coordinates
#   a: linear part of curve: y^2 = x^3 + ax + b
#   b: constant part of curve
#   G: a curve point (G.x,G.y) used as a "generator" (starting point)
class ECcurve:
	def __init__(self):
		return

	# Prime field multiplication: return j*k mod p
	def field_mul(self,j,k):
		return (j*k)%self.p

	# Prime field division: return num/den mod p
	def field_div(self,num,den):
		inverse_den=modular_inverse(den%self.p,self.p)
		return self.field_mul(num%self.p,inverse_den)

	# Prime field exponentiation: raise num to power mod p
	def field_exp(self,num,power):
		return pow(num%self.p,power,self.p)

	# Return the special identity point
	#   We pick x=p, y=0
	def identity(self):
		return ECpoint(self,self.p,0)

	# Return true if point Q lies on our curve
	def touches(self,Q):
		y2=self.field_exp(Q.y,2)
		x3ab=(self.field_mul((Q.x*Q.x)%self.p+self.a,Q.x)+self.b)%self.p
		return y2==x3ab

	# Return the slope of the tangent of this curve at point Q
	def tangent(self,Q):
		return self.field_div(Q.x*Q.x*3+self.a,Q.y*2)

	# Return the (x,y) point where this line intersects our curve
	#  Q1 and Q2 are two points on the line of slope m
	def line_intersect(self,Q1,Q2,m):
		v=(Q1.y + self.p - (m*Q1.x)%self.p)%self.p
		x=(m*m + self.p-Q1.x + self.p-Q2.x)%self.p
		y=(self.p-(m*x)%self.p + self.p-v)%self.p
		return ECpoint(self,x,y)

	# Return a doubled version of this elliptic curve point
	def double(self,Q):
		if (Q.x==self.p): # doubling the identity
			return Q
		return self.line_intersect(Q,Q,self.tangent(Q))

	# Return the "sum" of these elliptic curve points
	def add(self,Q1,Q2):
		# Identity special cases
		if (Q1.x==self.p): # Q1 is identity
			return Q2
		if (Q2.x==self.p): # Q2 is identity
			return Q1

		# Equality special cases
		if (Q1.x==Q2.x): 
			if (Q1.y==Q2.y): # adding point to itself
				return self.double(Q1)
			else: # vertical pair--result is the identity
				return self.identity()

		# Ordinary case
		m=self.field_div(Q1.y+self.p-Q2.y,Q1.x+self.p-Q2.x)
		return self.line_intersect(Q1,Q2,m)

	# "Multiply" this elliptic curve point Q by the integer m
	#    Often the point Q will be the generator G
	def mul(self,Q,m):
		R=self.identity() # return point
		while m!=0:  # binary multiply loop
			if m&1: # bit is set
				# print("  mul: adding Q to R =",R);
				R=self.add(R,Q)
			m=m>>1
			if (m!=0):
				# print("  mul: doubling Q =",Q);
				Q=self.double(Q)
		
		return R

To run this, we need an elliptic curve. y^2 = x^3 + 7 is used by BitCoin, and taken from the "Standards for Efficient Cryptography (SEC)" group, named SECp256k1.

# This is the BitCoin elliptic curve, SECG secp256k1
#   See http://www.secg.org/SEC2-Ver-1.0.pdf
secp256k1=ECcurve()

# All arithmetic is modulo this big prime number
secp256k1.p=hex2int("FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F");
secp256k1.a=0 # it's a Koblitz curve, with no linear part
secp256k1.b=7 

# The "generator" G is a starting point on the curve
secp256k1.G=ECpoint(curve=secp256k1,
  x = hex2int("79BE667E F9DCBBAC 55A06295 CE870B07 029BFCDB 2DCE28D9 59F2815B 16F81798"),
  y = hex2int("483ADA77 26A3C465 5DA4FBFC 0E1108A8 FD17B448 A6855419 9C47D08F FB10D4B8")
);

Finally, we might do something like Elliptic Curve Diffie-Hellman (ECDH) key exchange using this elliptic curve, like this:

# Diffie-Hellman key exchange
A_secret=rand.getrandbits(256)%curve.p
A_public=curve.mul(curve.G,A_secret);

B_secret=rand.getrandbits(256)%curve.p
B_public=curve.mul(curve.G,B_secret);

# A and B would exchange public points here

AB_shared=curve.mul(A_public,B_secret); # B computes this
BA_shared=curve.mul(B_public,A_secret); # A computes this

if (AB_shared.x == BA_shared.x):
	print("ECDH key exchange success! ",AB_shared);

Because elliptic curve point addition is commutative, multiplying the generator point by A's secret then B's, or by B's secret then A's, results in the same point.

Download my full Python3 elliptic curve point multiplication code here.