CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In
Avatar for 18.783 Fall 2023.

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download

Implementation of the group law on an Edwards curve

Views: 112
License: GPL3
Image: ubuntu2204
Kernel: SageMath 9.8

For the sake of comparison we first give formulas for Weierstrass curves in projective coordinates; scroll down to see the formulas for Edwards curves

def add_affine(P,Q,A):
    """ addition algorithm for affine points on an elliptic curve in short Weierstrass form y^2=x^3+Ax+B"""
    if P == O: return Q
    if Q == O: return P
    assert P[2] == 1 and Q[2] == 1  # we are using affine formulas, so make sure points are in affine form
    x1=P[0]; y1=P[1]; x2=Q[0]; y2=Q[1];
    if x1 != x2:
        m = (y2-y1)/(x2-x1)         # usual case: P and Q are distinct and not opposites
    else:
        if y1 == -y2: return O      # P = -Q (includes case where P=Q is 2-torsion)
        m = (3*x1^2+A) / (2*y1)     # P = Q so we are doubling
    x3 = m^2-x1-x2
    y3 = m*(x1-x3)-y1
    return (x3,y3,1)

def dbl_affine(P,A):
    """doubles a projective non-2-torsion point on an elliptic curve in short Weierstrass form"""
    x1=P[0]; y1=P[1];
    m = (3*x1^2+A) / (2*y1)
    x3 = m^2-2*x1
    y3 = m*(x1-x3)-y1
    return (x3,y3,1)

def add_projective(P,Q,A):
    """addition algorithm for projective points on an elliptic curve in short Weierstrass form y^2=x^3+Ax+B"""
    if P[2] == 0: return Q
    if Q[2] == 0: return P
    x1=P[0]; y1=P[1]; z1=P[2]; x2=Q[0]; y2=Q[1]; z2=Q[2]
    # We use formulas from http://hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#addition-add-1998-cmo-2
    x1z2 = x1*z2; x2z1 = x2*z1; y1z2 = y1*z2; y2z1 = y2*z1
    if x1z2 != x2z1:
        # usual case: P and Q are distinct and not opposites
        z1z2 = z1*z2
        u = y2z1-y1z2; uu = u^2
        v = x2z1-x1z2; vv = v^2; vvv = v^3
        r = vv*x1z2
        s = uu*z1z2 - vvv - 2*r
        x3 = v*s
        y3 = u*(r-s)-vvv*y1z2
        z3 = vvv*z1z2
    else:
        if y1z2 == -y2z1: return (0,1,0)      # P = -Q (includes case where P=Q is 2-torsion)
        # P = Q so we are doubling
        xx = x1^2; yy = y1^2; zz = z1^2
        w = A*zz + 3*xx
        s = 2*y1*z1; ss = s^2; sss = s*ss
        r = y1*s; rr = r^2
        b = (x1+r)^2-xx-rr
        h = w^2-2*b
        x3 = h*s
        y3 = w*(b-h)-2*rr
        z3 = sss
    return (x3,y3,z3)

def negate(P):
    """negates the point P on an elliptic curve in short Weierstrass form"""
    if P[2] == 0: return (0,1,0)
    return (P[0],-P[1],P[2])
    
def dbl_projective(P,A):
    """doubles a projective non-2-torsion point on an elliptic curve in short Weierstrass form"""
    x1 = P[0]; y1 = P[1]; z1 = P[2]
    xx = x1^2; yy = y1^2; zz = z1^2
    w = A*zz + 3*xx
    s = 2*y1*z1; ss = s^2; sss = s*ss
    r = y1*s; rr = r^2
    b = (x1+r)^2-xx-rr
    h = w^2-2*b
    x3 = h*s
    y3 = w*(b-h)-2*rr
    z3 = sss
    return (x3,y3,z3)

def equal_projective(P,Q):
    """Test whether two projective points are equal or not"""
    return P[0]*Q[2] == Q[0]*P[2] and P[1]*Q[2] == Q[1]*P[2]
    
def random_point(A,B):
    """generates a random affine point on the elliptic curve y^2 = x^2 + Ax + B"""
    F=A.parent()
    x0=F.random_element(); t=(x0^3+A*x0+B)
    while not t.is_square():
        x0=F.random_element(); t=(x0^3+A*x0+B)
    return (x0,sqrt(t),1)

p=random_prime(2^256,lbound=2^255)
F=GF(p)
A=F(3); B=F.random_element()
P0=random_point(A,B)
Q0=random_point(A,B)
print("Timing add_affine on random points...")
%timeit -n 30000 R=add_affine(P0,Q0,A)
print("Timing add_add_affine on equal points...")
%timeit -n 30000 R=add_affine(P0,P0,A)
print("Testing dbl_affine...")
%timeit -n 30000 R=dbl_affine(P0,A)
print("Timing add_projective on random points...")
%timeit -n 30000 R=add_projective(P0,Q0,A)
print("Timing add_projective on equal points...")
%timeit -n 30000 R=add_projective(P0,P0,A)
print("Testing dbl_projective...")
%timeit -n 30000 R=dbl_projective(P0,A)
Timing add_affine on random points... 6.34 µs ± 1.02 µs per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_add_affine on equal points... 8.06 µs ± 1.01 µs per loop (mean ± std. dev. of 7 runs, 30000 loops each) Testing dbl_affine... 6.97 µs ± 402 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_projective on random points... 11.5 µs ± 1.7 µs per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_projective on equal points... 16.7 µs ± 2.83 µs per loop (mean ± std. dev. of 7 runs, 30000 loops each) Testing dbl_projective... 11.2 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 30000 loops each)

Note that distinct triples of proejctive coordinates may define the same projective point. To compare points for equality we need to use the equal_projective function above (which also works for Edwards curves)

R1=add_projective(P0,Q0,A)
R2=add_projective(Q0,P0,A)
print("R1 = " + str(R1))
print("R2 = " + str(R2))
print("R1 == R2: " + str(R1 == R2))
print("equal_projective(R1,R2): %s" % equal_projective(R1,R2))
R1 = (39631903559207355134813242690485305706327172391023199616984319804673277252824, 33869231864040439687106036648793069871874765452794702250861206505319640701175, 1271603232244282005892045156571757933303549172868729422943298126929250743038) R2 = (39003209298511745779211371247669215406384172665513230424454802014897447644459, 44765880993678661226918577289361451240836579603741727790577915314251084196108, 77363509625474818908132568781582763179407795883667700618495823692641474154245) R1 == R2: False equal_projective(R1,R2): True 
def add_edwards(P,Q,d):
    """adds two projective points on an elliptic curve in Edwards form (assumes d non-square, so operation is complete)"""
    # uses unoptimized formula derived in class (could save 1M with optimization)
    x1=P[0]; y1=P[1]; z1=P[2]; x2=Q[0]; y2=Q[1]; z2=Q[2]
    x1y2 = x1*y2; x2y1 = x2*y1
    r = z1*z2; rr = r^2; s = x1y2+x2y1; t = d*x1y2*x2y1; u = y1*y2 - x1*x2
    v = rr+t; w = rr-t
    x3 = r*s*w
    y3 = r*u*v
    z3 = v*w
    return (x3,y3,z3)
    
def dbl_edwards(P):
    """doubles a projective point on an elliptic curve in Edwards form"""
    # uses optimized formula from Bernstein-Lange
    x1 = P[0]; y1 = P[1]; z1 = P[2]
    B=(x1+y1)^2; C=x1^2; D=y1^2; E = C+D;
    J = E-2*z1^2; x3 = (B-E)*J; y3 = E*(C-D); z3 = E*J
    return (x3,y3,z3)
    
def neg_edwards(P):
    return (-P[0],P[1],P[2])
    
def random_edwards_point(d):
    """generates a random projective point on the Edwards curve x^2 + y^2 = 1+d*x^2*y^2 (assumes d non-square)"""
    F=d.parent()
    x0=F.random_element(); t=(1-x0^2)/(1-d*x0^2)
    while not t.is_square():
        x0=F.random_element(); t=(1-x0^2)/(1-d*x0^2)
    return (x0,sqrt(t),1)
    
def on_edwards_curve (P,d):
    x0=P[0]; y0=P[1]; z0=P[2]
    return z0^2*(x0^2+y0^2) == z0^4 + d*x0^2*y0^2
    
def is_identity_edwards (P):
    if P[0] == 0 and P[1] == P[2]: return true
    else: return false

p = 17
F=GF(p)
d = F.multiplicative_generator()
print("Affine points on Edwards curve x^2 + y^2 = 1 + %dx^2y^2 over F_%d"%(d,p))
for a in F:
    for b in F:
        if on_edwards_curve([a,b,1],d): print([a,b])
Affine points on Edwards curve x^2 + y^2 = 1 + 3x^2y^2 over F_17 [0, 1] [0, 16] [1, 0] [2, 5] [2, 12] [3, 4] [3, 13] [4, 3] [4, 14] [5, 2] [5, 15] [7, 7] [7, 10] [10, 7] [10, 10] [12, 2] [12, 15] [13, 3] [13, 14] [14, 4] [14, 13] [15, 5] [15, 12] [16, 0] 
p=random_prime(2^256,lbound=2^255)
F=GF(p)
d=F(2);
while d.is_square(): d+=1;
P0=random_edwards_point(d)
Q0=random_edwards_point(d)
print("Timing add_edwards on random points...")
%timeit -n 30000 add_edwards(P0,Q0,d)
print("Timing add_edwards on equal points...")
%timeit -n 30000 add_edwards(P0,P0,d)
print("Timing dbl_edwards...")
%timeit -n 30000 dbl_edwards(P0)
Timing add_edwards on random points... 9.31 µs ± 2.12 µs per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_edwards on equal points... 7.75 µs ± 1.14 µs per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing dbl_edwards... 4.86 µs ± 541 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)