Edwards curves

An elliptic curve in Edwards form is a curve of the form x^2 + y^2 = c^2(1 + dx^2y^2). where cd(1-dc^4) is nonzero. (Technically, the curve is not elliptic, because it has singularities; but resolving those singularities produces an elliptic curve.) The neutral element of the curve is the point (0,c).

An Edwards-form elliptic curve x^2 + y^2 = c^2(1 + dx^2y^2) is birationally equivalent to the Montgomery-form elliptic curve (1/e)v^2 = u^3 + (4/e-2)u^2 + u where e = 1-dc^4. A point (x,y) on the Edwards curve, with nonzero x, corresponds to the point (u,v) on the Montgomery curve defined by u=(c+y)/(c-y) and v=2cu/x.

Edwards coordinates represent an affine point (x,y) on an Edwards-form elliptic curve x^2 + y^2 = c^2(1 + dx^2y^2) as (X:Y:Z) satisfying (X^2 + Y^2)Z^2 = c^2(Z^4 + dX^2Y^2). Here (X:Y:Z) = (sX:sY:sZ) for all nonzero s.

Speed records:

Edwards addition with c=1: 10M+1S+1D+7add. Algorithm: 2007 Bernstein/Lange.
Edwards strongly unified addition with c=1: 10M+1S+1D+7add. Algorithm: 2007 Bernstein/Lange.
Edwards mixed addition with c=1: 9M+1S+1D+7add. Algorithm: 2007 Bernstein/Lange.
Edwards addition with X2=1 with c=1: 9M+1S+1D+4add. Algorithm: 2007 Hisil/Carter/Dawson.
Edwards addition with Z1=1 and Z2=1 and c=1: 6M+1S+1D+8add. Algorithm: 2007 Bernstein/Lange.
Edwards doubling with c=1: 3M+4S+5add+1times2. Algorithm: 2007 Bernstein/Lange.
Edwards doubling with Z1=1: with c=1: 3M+3S+5add+1times2. Algorithm: 2007 Bernstein/Lange.
Edwards tripling with c=1: 9M+4S+6add+2times2. Algorithm: 2007 Bernstein/Birkner/Lange/Peters; independently 2007 Hisil/Carter/Dawson.

Some speed records are different when S/M is small:

Edwards addition with c=1: 7M+5S+1D+13add+2times2. Algorithm: 2007 Bernstein/Lange.
Edwards strongly unified addition with c=1: 7M+5S+1D+13add+2times2. Algorithm: 2007 Bernstein/Lange.
Edwards mixed addition with c=1: 6M+5S+1D+13add+2times2. Algorithm: 2007 Bernstein/Lange.
Edwards tripling with c=1: 7M+7S+12add+2times2+1times4. Algorithm: 2007 Bernstein/Birkner/Lange/Peters.

The following commands for the Magma computer-algebra system check various addition formulas for Edwards coordinates.

Edwards scaling.

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     // here are the formulas:
     A:=1/Z1;
     X2:=X1*A;
     Y2:=Y1*A;
     Z2:=1;
     // check:
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     S!(Z2-1);
     S!(x2-x1); S!(y2-y1);

Edwards addition (10M+1S+1C+1D+7add) matches traditional addition. 2007 Edwards, strongly unified:

     K<c,d,x1,x2>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<y1,y2>:=PolynomialRing(K,2);
     S:=quo<R|x1^2+y1^2-c^2*(1+d*x1^2*y1^2),x2^2+y2^2-c^2*(1+d*x2^2*y2^2)>;
     // the Edwards addition law:
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2));
     y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     // map to the Montgomery curve:
     u1:=(c+y1)/(c-y1); v1:=2*c*u1/x1; S!((1/e)*v1^2-u1^3-(4/e-2)*u1^2-u1);
     u2:=(c+y2)/(c-y2); v2:=2*c*u2/x2; S!((1/e)*v2^2-u2^3-(4/e-2)*u2^2-u2);
     u3:=(c+y3)/(c-y3); v3:=2*c*u3/x3; S!((1/e)*v3^2-u3^3-(4/e-2)*u3^2-u3);
     // add on the Montgomery curve:
     lambda:=(v2-v1)/(u2-u1);
     r3:=(1/e)*lambda^2-(4/e-2)-u1-u2; s3:=lambda*(u1-r3)-v1;
     // check the answer:
     S!(u3-r3); S!(v3-s3);

2007 Bernstein/Lange, 10M + 1S + 1C + 1D + 7add, strongly unified:

     K<c,d,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),6));
     e:=1-d*c^4;
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     A:=Z1*Z2;
     B:=A^2;
     C:=X1*X2;
     D:=Y1*Y2;
     E:=d*C*D;
     F:=B-E;
     G:=B+E;
     X3:=A*F*((X1+Y1)*(X2+Y2)-C-D);
     Y3:=A*G*(D-C);
     Z3:=c*F*G;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 10M + 1S + 1C + 1D + 7add, strongly unified, two temporary registers:

     K<c,d,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),6));
     e:=1-d*c^4;
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     R1:=X1; R2:=Y1; R3:=Z1;
     R4:=X2; R5:=Y2; R6:=Z2;
     R3:=R3*R6;
     R7:=R1+R2;
     R8:=R4+R5;
     R1:=R1*R4;
     R2:=R2*R5;
     R7:=R7*R8;
     R7:=R7-R1;
     R7:=R7-R2;
     R7:=R7*R3;
     R8:=R1*R2;
     R8:=d*R8;
     R2:=R2-R1;
     R2:=R2*R3;
     R3:=R3^2;
     R1:=R3-R8;
     R3:=R3+R8;
     R2:=R2*R3;
     R3:=R3*R1;
     R1:=R1*R7;
     R3:=c*R3;
     X3:=R1; Y3:=R2; Z3:=R3;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 7M + 5S + 1C + 1D + 13add + 2times2, strongly unified:

     K<c,d,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),6));
     e:=1-d*c^4;
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     A:=Z1*Z2;
     B:=A^2;
     C:=X1*X2;
     D:=Y1*Y2;
     E:=d*C*D;
     BB:=B^2;
     EE:=E^2;
     H:=(A+B)^2-BB;
     I:=(A+E)^2-EE;
     X3:=(H-I)*((X1+Y1)*(X2+Y2)-C-D);
     Y3:=(H+I-2*B)*(D-C);
     Z3:=2*c*(BB-EE);
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 10M + 1S + 1C + 1D + 9add + 3timesi + 2times2, strongly unified:

     Qpoly<i>:=PolynomialRing(Rationals());
     Gauss<i>:=quo<Qpoly|i^2+1>;
     K<c,d,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Gauss,6));
     e:=1-d*c^4;
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     iX2:=i*X2;
     C2:=Y2+iX2;
     D2:=Y2-iX2;
     iX1:=i*X1;
     C1:=Y1+iX1;
     D1:=Y1-iX1;
     A:=Z1*Z2;
     B:=2*A^2;
     C:=C1*C2;
     D:=D1*D2;
     L:=D+C;
     M:=Y1*Y2;
     N:=2*M-L;
     E:=d*M*N;
     F:=B-E;
     G:=B+E;
     X3:=i*A*F*(D-C);
     Y3:=A*G*L;
     Z3:=c*G*F;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

Edwards mixed addition (9M+1S+1C+1D+7add) matches traditional addition. 2007 Bernstein/Lange, 9M + 1S + 1C + 1D + 7add:

     K<c,d,X1,Y1,X2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     e:=1-d*c^4;
     R<Z1,Y2>:=PolynomialRing(K,2);
     Z2:=1;
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     B:=Z1^2;
     C:=X1*X2;
     D:=Y1*Y2;
     E:=d*C*D;
     F:=B-E;
     G:=B+E;
     X3:=Z1*F*((X1+Y1)*(X2+Y2)-C-D);
     Y3:=Z1*G*(D-C);
     Z3:=c*F*G;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 9M + 1S + 1C + 1D + 7add, two temporary registers:

     K<c,d,X1,Y1,X2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     e:=1-d*c^4;
     R<Z1,Y2>:=PolynomialRing(K,2);
     Z2:=1;
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     R1:=X1; R2:=Y1; R3:=Z1;
     R4:=X2; R5:=Y2;
     R7:=R1+R2;
     R6:=R4+R5;
     R1:=R1*R4;
     R2:=R2*R5;
     R7:=R7*R6;
     R7:=R7-R1;
     R7:=R7-R2;
     R7:=R7*R3;
     R6:=R1*R2;
     R6:=d*R6;
     R2:=R2-R1;
     R2:=R2*R3;
     R3:=R3^2;
     R1:=R3-R6;
     R3:=R3+R6;
     R2:=R2*R3;
     R3:=R3*R1;
     R1:=R1*R7;
     R3:=c*R3;
     X3:=R1; Y3:=R2; Z3:=R3;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 6M + 5S + 1C + 1D + 13add + 2times2:

     K<c,d,X1,Y1,X2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     e:=1-d*c^4;
     R<Z1,Y2>:=PolynomialRing(K,2);
     Z2:=1;
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     B:=Z1^2;
     C:=X1*X2;
     D:=Y1*Y2;
     E:=d*C*D;
     BB:=B^2;
     EE:=E^2;
     H:=(Z1+B)^2-BB;
     I:=(Z1+E)^2-EE;
     X3:=(H-I)*((X1+Y1)*(X2+Y2)-C-D);
     Y3:=(H+I-2*B)*(D-C);
     Z3:=2*c*(BB-EE);
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

Edwards addition with X2=1 (9M + 1S + 1C + 1D + 4add) matches traditional addition. 2007 Hisil/Carter/Dawson, 9M + 1S + 1C + 1D + 4add, two temporary registers:

     K<c,d,X1,Y1,Z2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     e:=1-d*c^4;
     R<Z1,Y2>:=PolynomialRing(K,2);
     X2:=1;
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     T0:=X1*Y2;
     T0:=T0+Y1;
     Y3:=Y1*Y2;
     T1:=Y3*X1;
     Y3:=Y3-X1;
     Z3:=Z1*Z2;
     X3:=T0*Z3;
     Y3:=Y3*Z3;
     T1:=d*T1;
     Z3:=Z3^2;
     T0:=Z3-T1;
     Z3:=Z3+T1;
     X3:=X3*T0;
     Y3:=Y3*Z3;
     Z3:=Z3*T0;
     Z3:=c*Z3;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

Edwards addition with Z1=1 and Z2=1 (6M+1S+1C+1D+8add) matches traditional addition. 2007 Bernstein/Lange, 6M + 1S + 1C + 1D + 8add:

     K<c,d,X1,X2>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Y1,Y2>:=PolynomialRing(K,2);
     Z1:=1;
     Z2:=1;
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2),(X2^2+Y2^2)*Z2^2-c^2*(Z2^4+d*X2^2*Y2^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     C:=X1*X2;
     D:=Y1*Y2;
     E:=d*C*D;
     X3:=(1-E)*((X1+Y1)*(X2+Y2)-C-D);
     Y3:=(1+E)*(D-C);
     Z3:=c*(1-E^2);
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

Edwards doubling (3M+4S+3C+5add+1times2) matches traditional doubling. 2007 Edwards, strongly unified:

     K<c,d,x1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     e:=1-d*c^4;
     R<y1>:=PolynomialRing(K,1);
     S:=quo<R|x1^2+y1^2-c^2*(1+d*x1^2*y1^2)>;
     x2:=x1; y2:=y1;
     // the Edwards addition law:
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2));
     y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     // map to the Montgomery curve:
     u1:=(c+y1)/(c-y1); v1:=2*c*u1/x1; S!((1/e)*v1^2-u1^3-(4/e-2)*u1^2-u1);
     u2:=(c+y2)/(c-y2); v2:=2*c*u2/x2; S!((1/e)*v2^2-u2^3-(4/e-2)*u2^2-u2);
     u3:=(c+y3)/(c-y3); v3:=2*c*u3/x3; S!((1/e)*v3^2-u3^3-(4/e-2)*u3^2-u3);
     // double on the Montgomery curve:
     lambda:=(3*u1^2+2*(4/e-2)*u1+1)/((2/e)*v1);
     r3:=(1/e)*lambda^2-(4/e-2)-u1-u2; s3:=lambda*(u1-r3)-v1;
     // check the answer:
     S!(u3-r3); S!(v3-s3);

2007 Bernstein/Lange, 10M + 1S + 1C + 1D + 7add, strongly unified:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=x1; y2:=y1;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // here are the formulas:
     A:=Z1*Z2;
     B:=A^2;
     C:=X1*X2;
     D:=Y1*Y2;
     E:=d*C*D;
     F:=B-E;
     G:=B+E;
     X3:=A*F*((X1+Y1)*(X2+Y2)-C-D);
     Y3:=A*G*(D-C);
     Z3:=c*F*G;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 10M + 1S + 1C + 1D + 7add, two temporary registers, strongly unified:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=x1; y2:=y1;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     R1:=X1; R2:=Y1; R3:=Z1;
     R4:=X1; R5:=Y1; R6:=Z1;
     R3:=R3*R6;
     R7:=R1+R2;
     R8:=R4+R5;
     R1:=R1*R4;
     R2:=R2*R5;
     R7:=R7*R8;
     R7:=R7-R1;
     R7:=R7-R2;
     R7:=R7*R3;
     R8:=R1*R2;
     R8:=d*R8;
     R2:=R2-R1;
     R2:=R2*R3;
     R3:=R3^2;
     R1:=R3-R8;
     R3:=R3+R8;
     R2:=R2*R3;
     R3:=R3*R1;
     R1:=R1*R7;
     R3:=c*R3;
     X3:=R1; Y3:=R2; Z3:=R3;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 7M + 5S + 1C + 1D + 13add + 2times2, strongly unified:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=x1; y2:=y1;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // here are the formulas:
     A:=Z1*Z2;
     B:=A^2;
     C:=X1*X2;
     D:=Y1*Y2;
     E:=d*C*D;
     BB:=B^2;
     EE:=E^2;
     H:=(A+B)^2-BB;
     I:=(A+E)^2-EE;
     X3:=(H-I)*((X1+Y1)*(X2+Y2)-C-D);
     Y3:=(H+I-2*B)*(D-C);
     Z3:=2*c*(BB-EE);
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 3M + 4S + 3C + 5add + 1times2:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=x1; y2:=y1;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     B:=(X1+Y1)^2;
     C:=X1^2;
     D:=Y1^2;
     E:=C+D;
     H:=(c*Z1)^2;
     J:=E-2*H;
     X3:=c*(B-E)*J;
     Y3:=c*E*(C-D);
     Z3:=E*J;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 3M + 4S + 3C + 5add + 1times2, two temporary registers:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=x1; y2:=y1;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     R1:=X1; R2:=Y1; R3:=Z1;
     R4:=R1+R2;
     R3:=c*R3;
     R1:=R1^2;
     R2:=R2^2;
     R3:=R3^2;
     R4:=R4^2;
     R3:=2*R3;
     R5:=R1+R2;
     R2:=R1-R2;
     R4:=R4-R5;
     R3:=R5-R3;
     R1:=R3*R4;
     R3:=R3*R5;
     R2:=R2*R5;
     R1:=c*R1;
     R2:=c*R2;
     X3:=R1; Y3:=R2; Z3:=R3;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Lange, 3M + 4S + 3C + 5add + 2times2, one temporary register:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=x1; y2:=y1;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     R1:=X1; R2:=Y1; R3:=Z1;
     R3:=c*R3;
     R4:=R1^2;
     R1:=R1+R2;
     R1:=R1^2;
     R2:=R2^2;
     R3:=R3^2;
     R3:=2*R3;
     R4:=R2+R4;
     R2:=2*R2;
     R2:=R4-R2;
     R1:=R1-R4;
     R2:=R2*R4;
     R3:=R4-R3;
     R1:=R1*R3;
     R3:=R3*R4;
     R1:=c*R1;
     R2:=c*R2;
     X3:=R1; Y3:=R2; Z3:=R3;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

Edwards doubling with Z1=1 (3M+3S+3C+5add+1times2) matches traditional doubling. 2007 Bernstein/Lange, 3M + 3S + 3C + 5add + 1times2:

     K<c,d,X1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     e:=1-d*c^4;
     Z1:=1;
     R<Y1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^2+y1^2-c^2*(1+d*x1^2*y1^2));
     x2:=x1; y2:=y1;
     S!(x2^2+y2^2-c^2*(1+d*x2^2*y2^2));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2)); y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     S!(x3^2+y3^2-c^2*(1+d*x3^2*y3^2));
     // here are the formulas:
     B:=(X1+Y1)^2;
     C:=X1^2;
     D:=Y1^2;
     E:=C+D;
     H:=c*c;
     J:=E-2*H;
     X3:=c*(B-E)*J;
     Y3:=c*E*(C-D);
     Z3:=E*J;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

Edwards tripling for c=1 (9M+4S+6add+2times2) matches traditional tripling. S,M counts stated in the literature: 9M + 4S in 2007 Bernstein/Birkner/Lange/Peters. 9M + 4S in 2007 Hisil/Carter/Dawson, independently of 2007 Bernstein/Birkner/Lange/Peters. 7M + 7S in 2007 Bernstein/Birkner/Lange/Peters.

2007 Hisil/Carter/Dawson, 9M + 4S + 1C + 13add + 2times2:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     x2:=(x1*y1+y1*x1)/(c*(1+d*x1*x1*y1*y1));
     y2:=(y1*y1-x1*x1)/(c*(1-d*x1*x1*y1*y1));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2));
     y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     // here are the formulas:
     A := X1^2;
     B := Y1^2;
     C := (2*c*Z1)^2;
     D := (A+B)^2;
     E := 2*(A+B)*(A-B);
     F := A*C;
     G := B*C;
     X3 := X1*(E-(D-G))*(D-G);
     Y3 := Y1*(E+(D-F))*(D-F);
     Z3 := Z1*(E-(D-G))*(E+(D-F));
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Birkner/Lange/Peters, 9M + 4S + 1C + 6add + 2times2:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     x2:=(x1*y1+y1*x1)/(c*(1+d*x1*x1*y1*y1));
     y2:=(y1*y1-x1*x1)/(c*(1-d*x1*x1*y1*y1));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2));
     y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     // here are the formulas:
     XX:=X1^2;
     YY:=Y1^2;
     ZZ:=(2*c*Z1)^2;
     D:=XX+YY;
     DD:=D^2;
     H:=2*D*(XX-YY);
     P:=DD-YY*ZZ;
     Q:=DD-XX*ZZ;
     T:=H+Q;
     U:=H-P;
     X3:=P*U*X1;
     Y3:=Q*T*Y1;
     Z3:=T*U*Z1;
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Birkner/Lange/Peters, 7M + 7S + 12add + 2times2 + 1times4, if c=1:

     K<d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     c:=1;
     e:=1-d*c^4;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     x2:=(x1*y1+y1*x1)/(c*(1+d*x1*x1*y1*y1));
     y2:=(y1*y1-x1*x1)/(c*(1-d*x1*x1*y1*y1));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2));
     y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     // here are the formulas:
     XX:=X1^2;
     YY:=Y1^2;
     ZZ:=Z1^2;
     ZZ4:=4*ZZ;
     D:=XX+YY;
     DD:=D^2;
     H:=2*D*(XX-YY);
     P:=DD-YY*ZZ4;
     Q:=DD-XX*ZZ4;
     T:=H+Q;
     TT:=T^2;
     U:=H-P;
     X3:=2*P*U*X1;
     Y3:=Q*((T+Y1)^2-TT-YY);
     Z3:=U*((T+Z1)^2-TT-ZZ);
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2007 Bernstein/Birkner/Lange/Peters, 7M + 7S + 1C^2 + 12add + 2times2 + 1times4:

     K<c,d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     e:=1-d*c^4;
     a:=c^2;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|(X1^2+Y1^2)*Z1^2-c^2*(Z1^4+d*X1^2*Y1^2)>;
     x1:=X1/Z1; y1:=Y1/Z1;
     x2:=(x1*y1+y1*x1)/(c*(1+d*x1*x1*y1*y1));
     y2:=(y1*y1-x1*x1)/(c*(1-d*x1*x1*y1*y1));
     x3:=(x1*y2+y1*x2)/(c*(1+d*x1*x2*y1*y2));
     y3:=(y1*y2-x1*x2)/(c*(1-d*x1*x2*y1*y2));
     // here are the formulas:
     XX:=X1^2;
     YY:=Y1^2;
     ZZ:=Z1^2;
     ZZ4:=4*a*ZZ;
     D:=XX+YY;
     DD:=D^2;
     H:=2*D*(XX-YY);
     P:=DD-YY*ZZ4;
     Q:=DD-XX*ZZ4;
     T:=H+Q;
     TT:=T^2;
     U:=H-P;
     X3:=2*P*U*X1;
     Y3:=Q*((T+Y1)^2-TT-YY);
     Z3:=U*((T+Z1)^2-TT-ZZ);
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3);