I am reading the section about elliptic curve cryptography in Christoffer Paares book ("Understanding Cryptography"). I decided that I would implement a function for elliptic curve point addition and point doubling in python. For my test I used the example in the book so I could test my results.
The curve used in the example is: y^2 = x^3 + 2x + 2 mod 17
The generator used is: p = (5,1)
Thus the cycle becomes:
1p = (5,1)
2p = (6,3)
3p = (10,6)
4p = (3,1)
5p = (9,16)
6p = (16,13)
7p = (0,6)
8p = (13,7)
9p = (7,6)
10p = (7,1)
11p = (13,10)
12p = (0,11)
13p = (16,4)
14p = (9,1)
15p = (3,16)
16p = (10,11)
17p = (6,14)
18p = (5,16)
19p = The neutral element (Point at infinity)
20p = (5,1)
...
I wrote this code in an attempt to reproduce the results:
def add(a,p,P,Q):
#Check For Neutral Element
if P == (0,0) or Q == (0,0):
return (P[0]+Q[0],P[1]+Q[1])
#Check For Inverses (Return Neutral Element - Point At Infinity)
if P[0] == Q[0]:
if (-P[1])%p == Q[1] or (-Q[1])%p == P[1]:
return (0,0)
#Calculate Slope
if P != Q:
s = (Q[1]-P[1]) / (Q[0]-P[0])
else:
s = ((3*(P[0]*P[0])+a)%p) ** (2*P[1])
#Calculate Third Intersection
x = s*s - P[0] - Q[0]
y = (s*(P[0]-x)) - P[1]
y = y%p
x = x%p
return (x,y)
r = (5,1)
for i in range(1,20):
print i,':',r
r = add(2,17,r,(5,1))
However the output is:
As you might see it follows the expected result until 6p and then enters a cycle with the length of 2. I have been staring at this for hours now and I still don't know why it doesn't work (after all: how hard can it be... it's 30 lines of python) .
I am not really aware of the topic but here is a link to a repository implementing the ECC: github
Edit: The actual problem is the division modulo p. You can't just divide using integer arithmetics (15 / 4 == 3) but need to multiply by the inverse of 4 modulo 17 instead. The inverse of 4 modulo 17 is 13, because 4 * 13 % 17 == 1. Your fraction becomes 15*13, which is the equivalent of saying »15 * 1/4 modulo 17«. Just put some debug prints around your slope computation and you will see when the inversion begins to differ from simple integer division.
def euclid(a, b):
'''Solve x*a + y*b = ggt(a, b) and return (x, y, ggt(a, b))'''
# Non-recursive approach hence suitable for large numbers
x = yy = 0
y = xx = 1
while b:
q = a // b
a, b = b, a % b
x, xx = xx - q * x, x
y, yy = yy - q * y, y
return xx, yy, a
def inv(a, n):
'''Perform inversion 1/a modulo n. a and n should be COPRIME.'''
# coprimality is not checked here in favour of performance
i = euclid(a, n)[0]
while i < 0:
i += n
return i
def add(a,p,P,Q):
#Check For Neutral Element
if P == (0,0) or Q == (0,0):
return (P[0]+Q[0],P[1]+Q[1])
#Check For Inverses (Return Neutral Element - Point At Infinity)
if P[0] == Q[0]:
if (-P[1])%p == Q[1] or (-Q[1])%p == P[1]:
return (0,0)
#Calculate Slope
if P != Q:
# perform multiplication with the inverse modulo p
s = (Q[1]-P[1]) * inv(Q[0]-P[0], p)
else:
s = ((3*(P[0]*P[0])+a)%p) ** (2*P[1])
#Calculate Third Intersection
x = s*s - P[0] - Q[0]
y = (s*(P[0]-x)) - P[1]
y = y%p
x = x%p
return (x,y)
r = (5,1)
for i in range(1,20):
print i,':',r
r = add(2,17,r,(5,1))
prints
1 : (5, 1)
2 : (6, 3)
3 : (10, 6)
4 : (3, 1)
5 : (9, 16)
6 : (16, 13)
7 : (0, 6)
8 : (13, 7)
9 : (7, 6)
10 : (7, 11)
11 : (13, 10)
12 : (0, 11)
13 : (16, 4)
14 : (9, 1)
15 : (3, 16)
16 : (10, 11)
17 : (6, 14)
18 : (5, 16)
19 : (0, 0)
Here a link to pypi. Feel free to use or improve it.