Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

Explicit formulas for the group law on Edwards curves with comparison to formulas for Weierstrass curves in projective coordinates, as presented in 18.783 Lecture 2

770 views
License: AGPL3
unlisted
ubuntu2004
Tweet
Share
Share
Kernel: SageMath 9.2

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... 4.82 µs ± 344 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_add_affine on equal points... 6.73 µs ± 321 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Testing dbl_affine... 6.33 µs ± 280 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_projective on random points... 8.75 µs ± 178 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_projective on equal points... 12.6 µs ± 159 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Testing dbl_projective... 9.13 µs ± 65.7 ns 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... 6.44 µs ± 366 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing add_edwards on equal points... 6.46 µs ± 296 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each) Timing dbl_edwards... 4.96 µs ± 559 ns per loop (mean ± std. dev. of 7 runs, 30000 loops each)