Jacobi quartics

An elliptic curve in Jacobi-quartic form is a curve of the form y^2 = x^4 + 2a x^2 + 1 where a^2-1 is nonzero. The neutral element of the curve is the point (0,1).

A Jacobi-quartic elliptic curve y^2 = x^4 + 2a x^2 + 1 is birationally equivalent to the Weierstrass-form elliptic curve 2v^2 = u^3 - 2au^2 + (a^2-1)u. A typical point (x,y) on the Jacobi quartic corresponds to the point (u,v) on the Weierstrass curve defined by u = a + (y+1)/x^2, v = u/x.

Jacobi-quartic coordinates represent an affine point (x,y) on a Jacobi quartic y^2 = x^4 + 2a x^2 + 1 as (X:Y:Z) satisfying Y^2 = X^4 + 2a X^2 Z^2 + Z^4. Here (X:Y:Z) = (sX:s^2 Y:sZ) for all nonzero s.

Speed records:

Jacobi-quartic addition: 10M+3S+1D+13add+1times2. Algorithm: 2001 Billet Joye.
Jacobi-quartic strongly unified addition: 10M+3S+1D+13add+1times2. Algorithm: 2001 Billet Joye.
Cached Jacobi-quartic strongly unified addition: 9M+3S+1D+12add+1times2 after 1M+1add for caching X2*Z2+Y2. Algorithm: 2001 Billet Joye.
Cached Jacobi-quartic strongly unified addition, alternative approach: 8M+3S+1D+11add+2times2+1times4 after 3S+4add for caching X2^2,2*X2*Z2,Z2^2,X2^2+Z2^2,X2^2+Z2^2+Y2. Algorithm: 2007 Duquesne, plus two S-M tradeoffs.
Jacobi-quartic mixed addition: 8M+3S+1D+13add+1times2. Algorithm: 2001 Billet Joye.
Jacobi-quartic addition with Z1=1 and Z2=1: 5M + 2S + 1D + 10add + 1times2. Algorithm: 2001 Billet Joye.
Jacobi-quartic doubling: 2M + 6S + 1D + 9add + 2times2. Algorithm: 2007 Feng/Wu. Alternative using alpha = sqrt(1-delta^2): 1M + 7S + 3D + 7add + 3times2. Algorithm: 2007 Feng/Wu.
Jacobi-quartic doubling with Z1=1: 1M + 4S + 1D + 9add + 1times2 + 1times4. Algorithm: 2007 Feng/Wu.

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

Jacobi-quartic scaling.

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

Jacobi-quartic negation.

     K<a,x1>:=FieldOfFractions(PolynomialRing(Rationals(),2));
     R<y1>:=PolynomialRing(K,1);
     S:=quo<R|y1^2-x1^4-2*a*x1^2-1>;
     // here are the formulas:
     x2:=-x1;
     y2:=y1;
     // map to the Weierstrass curve:
     u1:=a+(y1+1)/x1^2; v1:=u1/x1; S!(2*v1^2-u1^3+2*a*u1^2-(a^2-1)*u1);
     u2:=a+(y2+1)/x2^2; v2:=u2/x2; S!(2*v2^2-u2^3+2*a*u2^2-(a^2-1)*u2);
     // check the answer:
     S!(u2-u1); S!(v2+v1);

Jacobi-quartic addition (10M+3S+1D+13add+1times2) matches traditional addition. 2001 Billet Joye, formula (11), specialized to affine, strongly unified:

     K<a,x1,x2>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     R<y1,y2>:=PolynomialRing(K,2);
     S:=quo<R|y1^2-x1^4-2*a*x1^2-1,y2^2-x2^4-2*a*x2^2-1>;
     // here are the formulas:
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     // map to the Weierstrass curve:
     u1:=a+(y1+1)/x1^2; v1:=u1/x1; S!(2*v1^2-u1^3+2*a*u1^2-(a^2-1)*u1);
     u2:=a+(y2+1)/x2^2; v2:=u2/x2; S!(2*v2^2-u2^3+2*a*u2^2-(a^2-1)*u2);
     u3:=a+(y3+1)/x3^2; v3:=u3/x3; S!(2*v3^2-u3^3+2*a*u3^2-(a^2-1)*u3);
     // add on the Weierstrass curve:
     lambda:=(v2-v1)/(u2-u1);
     r3:=2*lambda^2+2*a-u1-u2; s3:=lambda*(u1-r3)-v1;
     // check the answer:
     S!(u3-r3); S!(v3-s3);

2001 Billet Joye, formula (11), strongly unified:

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

2001 Billet Joye, page 6, 10M + 3S + 1D + 13add + 1times2, 8 registers, strongly unified:

     K<a,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,Y2^2-X2^4-2*a*X2^2*Z2^2-Z2^4>;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     T1:=X1; T2:=Y1; T3:=Z1;
     T4:=X2; T5:=Y2; T6:=Z2;
     T7:=T1*T3;
     T7:=T2+T7;
     T8:=T4*T6;
     T8:=T5+T8;
     T2:=T2*T5;
     T7:=T7*T8;
     T7:=T7-T2;
     T5:=T1*T4;
     T1:=T1+T3;
     T8:=T3*T6;
     T4:=T4+T6;
     T6:=T5*T8;
     T7:=T7-T6;
     T1:=T1*T4;
     T1:=T1-T5;
     T1:=T1-T8;
     T3:=T1^2;
     T6:=2*T6;
     T3:=T3-T6;
     T3:=T3*T6;
     T4:=a*T6;
     T2:=T2+T4;
     T4:=T8^2;
     T8:=T5^2;
     T5:=T4+T8;
     T2:=T2*T5;
     T2:=T2+T3;
     T5:=T4-T8;
     X3:=T7; Y3:=T2; Z3:=T5;
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

9M + 3S + 1D + 12add + 1times2 after caching 1M + 1add, strongly unified:

     K<a,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,Y2^2-X2^4-2*a*X2^2*Z2^2-Z2^4>;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     // cached:
     X2Z2Y2:=X2*Z2+Y2;
     // uncached:
     T1:=X1; T2:=Y1; T3:=Z1;
     T4:=X2; T5:=Y2; T6:=Z2;
     T7:=T1*T3;
     T7:=T2+T7;
     T8:=X2Z2Y2;
     T2:=T2*T5;
     T7:=T7*T8;
     T7:=T7-T2;
     T5:=T1*T4;
     T1:=T1+T3;
     T8:=T3*T6;
     T4:=T4+T6;
     T6:=T5*T8;
     T7:=T7-T6;
     T1:=T1*T4;
     T1:=T1-T5;
     T1:=T1-T8;
     T3:=T1^2;
     T6:=2*T6;
     T3:=T3-T6;
     T3:=T3*T6;
     T4:=a*T6;
     T2:=T2+T4;
     T4:=T8^2;
     T8:=T5^2;
     T5:=T4+T8;
     T2:=T2*T5;
     T2:=T2+T3;
     T5:=T4-T8;
     X3:=T7; Y3:=T2; Z3:=T5;
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

2001 Billet Joye, formula (11), separating out X2^2,Z2^2 as suggested by 2007 Duquesne, 9M + 2S + 1D + 11add + 1times2 after 1M + 2S for caching X2^2, X2*Z2, Z2^2, strongly unified:

     K<a,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,Y2^2-X2^4-2*a*X2^2*Z2^2-Z2^4>;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     // cached:
     A2:=X2^2;
     B2:=X2*Z2;
     C2:=Z2^2;
     // uncached:
     A1:=X1^2;
     B1:=X1*Z1;
     C1:=Z1^2;
     A1A2:=A1*A2;
     B1B2:=B1*B2;
     C1C2:=C1*C2;
     Y1Y2:=Y1*Y2;
     E:=C1C2+A1A2;
     F:=(A1+C1)*(A2+C2)-E;
     G:=2*B1B2;
     X3:=(B1+Y1)*(B2+Y2)-B1B2-Y1Y2;
     Y3:=E*(Y1Y2+a*G)+G*F;
     Z3:=C1C2-A1A2;
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

2007 Duquesne, page 103, Table 1, 10M + 4S + 1D + 12add + 1times2:

     K<a,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     b:=-2*a;
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,Y2^2-X2^4-2*a*X2^2*Z2^2-Z2^4>;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     T1:=X1^2;
     T2:=X2^2;
     T3:=X1*Z1;
     T4:=X2*Z2;
     T5:=Z1^2;
     T6:=Z2^2;
     T7:=Y1;
     T8:=Y2;
     T9:=T7*T8;
     T7:=T7+T3;
     T8:=T8+T4;
     T3:=T3*T4;
     T7:=T7*T8;
     T7:=T7-T9;
     T7:=T7-T3;
     T4:=T1*T2;
     T8:=T5*T6;
     T1:=T1+T5;
     T2:=T2+T6;
     T5:=T1*T2;
     T5:=T5-T4;
     T5:=T5-T8;
     T1:=T8-T4;
     T2:=T8+T4;
     T6:=b*T3;
     T6:=T9-T6;
     T6:=T6*T2;
     T3:=2*T3;
     T3:=T5*T3;
     T8:=T6+T3;
     X3:=T7;
     Y3:=T8;
     Z3:=T1;
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

2007 Bernstein/Lange, 8M + 3S + 1D + 11add + 2times2 + 1times4 after 3S + 4add for caching X2^2, 2*X2*Z2, Z2^2, X2^2+Z2^2, X2^2+Z2^2+Y2, strongly unified:

     K<a,X1,Y1,X2,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,Y2^2-X2^4-2*a*X2^2*Z2^2-Z2^4>;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     // cached:
     A2:=X2^2;
     C2:=Z2^2;
     D2:=A2+C2;
     B2:=(X2+Z2)^2-D2;
     E2:=B2+Y2;
     // uncached:
     A1:=X1^2;
     C1:=Z1^2;
     D1:=A1+C1;
     B1:=(X1+Z1)^2-D1;
     E1:=B1+Y1;
     A1A2:=A1*A2;
     B1B2:=B1*B2;
     C1C2:=C1*C2;
     Y1Y2:=Y1*Y2;
     F:=C1C2+A1A2;
     G:=2*B1B2;
     X3:=E1*E2-B1B2-Y1Y2;
     Y3:=F*(4*Y1Y2+a*G)+(D1*D2-F)*G;
     Z3:=2*(C1C2-A1A2);
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

Jacobi-quartic mixed addition (8M+3S+1D+13add+1times2) matches traditional addition. 2001 Billet Joye, page 6, specialized to Z2=1, 8M + 3S + 1D + 13add + 1times2:

     K<a,X1,Y1,X2>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     R<Z1,Y2>:=PolynomialRing(K,2);
     Z2:=1;
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,Y2^2-X2^4-2*a*X2^2*Z2^2-Z2^4>;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     T1:=X1; T2:=Y1; T3:=Z1;
     T4:=X2; T5:=Y2;
     T7:=T1*T3;
     T7:=T2+T7;
     T8:=T4;
     T8:=T5+T8;
     T2:=T2*T5;
     T7:=T7*T8;
     T7:=T7-T2;
     T5:=T1*T4;
     T1:=T1+T3;
     T8:=T3;
     T4:=T4+1;
     T6:=T5*T8;
     T7:=T7-T6;
     T1:=T1*T4;
     T1:=T1-T5;
     T1:=T1-T8;
     T3:=T1^2;
     T6:=2*T6;
     T3:=T3-T6;
     T3:=T3*T6;
     T4:=a*T6;
     T2:=T2+T4;
     T4:=T8^2;
     T8:=T5^2;
     T5:=T4+T8;
     T2:=T2*T5;
     T2:=T2+T3;
     T5:=T4-T8;
     X3:=T7; Y3:=T2; Z3:=T5;
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

Jacobi-quartic addition with Z1=1 and Z2=1. 2001 Billet Joye, page 6, specialized to Z2=1, 5M + 2S + 1D + 10add + 1times2:

     K<a,X1,X2>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     R<Y1,Y2>:=PolynomialRing(K,2);
     Z1:=1; Z2:=1;
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,Y2^2-X2^4-2*a*X2^2*Z2^2-Z2^4>;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     T7:=Y1+X1;
     T8:=Y2+X2;
     T2:=Y1*Y2;
     T7:=T7*T8; 
     T7:=T7-T2; 
     T5:=X1*X2;
     T6:=T5;
     X3:=T7-T6;
     T1:=X1+X2; 
     T3:=T1^2; 
     T6:=2*T6; 
     T3:=T3-T6;
     T3:=T3*T6;
     T4:=a*T6; 
     T2:=T2+T4;
     T8:=T5^2; 
     T5:=T8+1; 
     T2:=T2*T5;
     Y3:=T2+T3;
     Z3:=1-T8;
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

Jacobi-quartic doubling (2M + 6S + 1D + 9add + 2times2) matches traditional doubling. 2001 Billet Joye, formula (11), specialized to affine, strongly unified:

     K<a,x1>:=FieldOfFractions(PolynomialRing(Rationals(),2));
     R<y1>:=PolynomialRing(K,1);
     S:=quo<R|y1^2-x1^4-2*a*x1^2-1>;
     x2:=x1; y2:=y1;
     // here are the formulas:
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     // map to the Weierstrass curve:
     u1:=a+(y1+1)/x1^2; v1:=u1/x1; S!(2*v1^2-u1^3+2*a*u1^2-(a^2-1)*u1);
     u2:=a+(y2+1)/x2^2; v2:=u2/x2; S!(2*v2^2-u2^3+2*a*u2^2-(a^2-1)*u2);
     u3:=a+(y3+1)/x3^2; v3:=u3/x3; S!(2*v3^2-u3^3+2*a*u3^2-(a^2-1)*u3);
     // add on the Weierstrass curve:
     lambda:=(3*u1^2-4*a*u1+(a^2-1))/(4*v1);
     r3:=2*lambda^2+2*a-u1-u2; s3:=lambda*(u1-r3)-v1;
     // check the answer:
     S!(u3-r3); S!(v3-s3);

2001 Billet Joye, formula (11), strongly unified:

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

2001 Billet Joye, formula (11), specialized to doubling:

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

2007 Bernstein/Lange, 1M + 9S + 1D + 10add + 2times2 + 1times4:

     K<a,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4>;
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     XX:=X1^2;
     XXXX:=XX^2;
     YY:=Y1^2;
     ZZ:=Z1^2;
     ZZZZ:=ZZ^2;
     M:=XX+ZZ;
     XZ:=(X1+Z1)^2-M;
     XZXZ:=XZ^2;
     X3:=(Y1+XZ)^2-YY-XZXZ;
     Y3:=(ZZZZ+XXXX)*(4*YY+2*a*XZXZ)+XZXZ^2;
     Z3:=2*(ZZZZ-XXXX);
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

2007 Hisil/Carter/Dawson, 2M + 6S + 2D + 5add + 1times2:

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

2007 Feng/Wu, first JQN2 doubling formula, obvious typos corrected:

     K<a,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     b:=4-4*a^2;
     epsilon:=1;
     delta:=-a;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4>;
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     U1:=(X1*Z1+Y1)^2;
     U2:=Y1^2;
     V1:=(X1^2)^2;
     S1:=(X1*Z1)^2;
     T:=U2-epsilon*V1+2*delta*S1;
     X3:=U1-U2-S1;
     Y3:=(T+epsilon*V1)*(U2-2*delta*S1)+4*epsilon*S1^2;
     Z3:=T-epsilon*V1;
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

2007 Feng/Wu, first JQN2 doubling formula, obvious typos corrected, epsilon=1, common subexpressions eliminated, 2M + 6S + 1D + 9add + 2times2:

     K<a,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     b:=4-4*a^2;
     delta:=-a;
     R<Z1>:=PolynomialRing(K,1);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4>;
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     XZ1:=X1*Z1;
     U1:=(XZ1+Y1)^2;
     U2:=Y1^2;
     V1:=(X1^2)^2;
     S1:=XZ1^2;
     deltaS1:=2*delta*S1;
     T:=U2-V1+deltaS1;
     X3:=U1-U2-S1;
     Y3:=(T+V1)*(U2-deltaS1)+(2*S1)^2;
     Z3:=T-V1;
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

2007 Feng/Wu, second JQN2 doubling formula, using alpha = sqrt(1 - delta^2):

     K<a,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     b:=4-4*a^2;
     delta:=-a;
     epsilon:=1;
     R<Z1,alpha>:=PolynomialRing(K,2);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,alpha^2+delta^2-epsilon>;
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     T1:=(X1*Z1)^2;
     T2:=Y1^2;
     S1:=(T2+2*alpha*T1)^2;
     U:=(X1*Z1+Y1)^2-T1-T2;
     V:=U^2;
     X3:=U;
     Y3:=S1-alpha*V;
     Z3:=T2+2*delta*T1-2*epsilon*X1^4;
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

2007 Feng/Wu, second JQN2 doubling formula, epsilon=1, using alpha = sqrt(1 - delta^2), common subexpressions removed, 1M + 7S + 3D + 7add + 3times2:

     K<a,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     b:=4-4*a^2;
     delta:=-a;
     epsilon:=1;
     R<Z1,alpha>:=PolynomialRing(K,2);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4,alpha^2+delta^2-epsilon>;
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     XX1:=X1^2;
     XZ1:=X1*Z1;
     T1:=XZ1^2;
     T2:=Y1^2;
     S1:=(T2+2*alpha*T1)^2;
     X3:=(XZ1+Y1)^2-T1-T2;
     V:=X3^2;
     Y3:=S1-alpha*V;
     Z3:=T2+2*delta*T1-2*XX1^2;
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);

Jacobi-quartic doubling with Z1=1. 2007 Feng/Wu, first JQN2 doubling formula, obvious typos corrected, epsilon=1, common subexpressions eliminated, 1M + 4S + 1D + 9add + 1times2 + 1times4:

     K<a,X1>:=FieldOfFractions(PolynomialRing(Rationals(),2));
     b:=4-4*a^2;
     delta:=-a; Z1:=1;
     R<Y1>:=PolynomialRing(K,1);
     S:=quo<R|Y1^2-X1^4-2*a*X1^2*Z1^2-Z1^4>;
     X2:=X1; Y2:=Y1; Z2:=Z1;
     // affine:
     x1:=X1/Z1; y1:=Y1/Z1^2; x2:=X2/Z2; y2:=Y2/Z2^2;
     r:=1-(x1*x2)^2;
     x3:=(x1*y2+y1*x2)/r;
     y3:=((1+(x1*x2)^2)*(y1*y2+2*a*x1*x2)+2*x1*x2*(x1^2+x2^2))/r^2;
     S!(y1^2-x1^4-2*a*x1^2-1);
     S!(y2^2-x2^4-2*a*x2^2-1);
     S!(y3^2-x3^4-2*a*x3^2-1);
     // here are the formulas:
     U1:=(X1+Y1)^2;
     U2:=Y1^2;
     S1:=X1^2;
     V1:=S1^2;
     deltaS1:=2*delta*S1;
     T:=U2-V1+deltaS1;
     X3:=U1-U2-S1;
     Y3:=(T+V1)*(U2-deltaS1)+4*V1;
     Z3:=T-V1;
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3^2);