Hessian curves

An elliptic curve in Hessian form is a curve of the form x^3 + y^3 + 1 = 3dxy where d^3-1 is nonzero. The neutral element of the curve is the point at infinity.

A Hessian-form elliptic curve x^3 + y^3 + 1 = 3dxy is birationally equivalent to the Weierstrass-form elliptic curve v^2 = u^3 - 27d(d^3+8)u + 54(d^6-20d^3-8). A typical point (x,y) on the Hessian curve corresponds to the point (u,v) on the Weierstrass curve defined by u = p - 9d^2 and v = 3p(y-x) where p = 12(d^3-1)/(d+x+y). This equivalence is based on the equivalence stated in 2001 Joye Quisquater, but is modified to match the standard choice of neutral element on the Hessian curve.

Hessian coordinates represent an affine point (x,y) on a Hessian-form elliptic curve x^3 + y^3 + 1 = 3dxy as (X:Y:Z) satisfying X^3 + Y^3 + Z^3 = 3dXYZ. Here (X:Y:Z) = (sX:sY:sZ) for all nonzero s.

Speed records:

Hessian addition: 12M + 3add. Algorithm: 1800s Sylvester.
Hessian weakly unified addition: 12M + 3add. Algorithm: 2001 Joye Quisquater.
Hessian readdition with X2=1: 5M + 6S + 12add after 3S + 3add + 2times2 to cache Y2^2, Z2^2, 2*Y2*Z2, 2*Y2, 2*Z2. Algorithm: 2007 Hisil/Carter/Dawson.
Hessian mixed addition: 10M + 3add. Algorithm: 1800s Sylvester.
Hessian addition with Z1=1 and Z2=1: 8M + 3add. Algorithm: 1800s Sylvester.
Hessian doubling: 7M + 1S + 8add. Algorithm: 2007 Hisil/Carter/Dawson.
Hessian doubling with Z1=1: 3M + 3S + 11add + 3times2. Algorithm: 2007 Bernstein/Lange.
Hessian tripling: 8M + 6S + 1D + 12add + 2times2. Algorithm: 2007 Hisil/Carter/Dawson.

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

Hessian doubling: 3M + 6S + 15add + 3times2. Algorithm: 2007 Hisil/Carter/Dawson.

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

Hessian scaling.

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

Hessian negation.

     K<d,x1>:=FieldOfFractions(PolynomialRing(Rationals(),2));
     R<y1>:=PolynomialRing(K,1);
     S:=quo<R|x1^3+y1^3+1-3*d*x1*y1>;
     // the Cauchy-Sylvester addition law:
     x2:=y1;
     y2:=x1;
     // map to the Weierstrass curve:
     p1:=12*(d^3-1)/(d+x1+y1); u1:=p1-9*d^2; v1:=3*p1*(y1-x1);
     S!(v1^2-u1^3+27*d*(d^3+8)*u1-54*(d^6-20*d^3-8));
     p2:=12*(d^3-1)/(d+x2+y2); u2:=p2-9*d^2; v2:=3*p2*(y2-x2);
     S!(v2^2-u2^3+27*d*(d^3+8)*u2-54*(d^6-20*d^3-8));
     // check the answer:
     S!(u2-u1); S!(v2+v1);

Hessian addition.
Addition law credited by 1986 Chudnovsky/Chudnovsky to Cauchy and Sylvester (see text before formula 4.21i):

     K<d,x1,x2>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     R<y1,y2>:=PolynomialRing(K,2);
     S:=quo<R|x1^3+y1^3+1-3*d*x1*y1,x2^3+y2^3+1-3*d*x2*y2>;
     // the Cauchy-Sylvester addition law:
     x3:=(y1^2*x2-y2^2*x1)/(x2*y2-x1*y1);
     y3:=(x1^2*y2-x2^2*y1)/(x2*y2-x1*y1);
     // map to the Weierstrass curve:
     p1:=12*(d^3-1)/(d+x1+y1); u1:=p1-9*d^2; v1:=3*p1*(y1-x1);
     S!(v1^2-u1^3+27*d*(d^3+8)*u1-54*(d^6-20*d^3-8));
     p2:=12*(d^3-1)/(d+x2+y2); u2:=p2-9*d^2; v2:=3*p2*(y2-x2);
     S!(v2^2-u2^3+27*d*(d^3+8)*u2-54*(d^6-20*d^3-8));
     p3:=12*(d^3-1)/(d+x3+y3); u3:=p3-9*d^2; v3:=3*p3*(y3-x3);
     S!(v3^2-u3^3+27*d*(d^3+8)*u3-54*(d^6-20*d^3-8));
     // add on the Weierstrass curve:
     lambda:=(v2-v1)/(u2-u1);
     r3:=lambda^2-u1-u2; s3:=lambda*(u1-r3)-v1;
     // check the answer:
     S!(u3-r3); S!(v3-s3);

Inversion-free version:

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

12M + 3add, consistent with the "12 multiplications" stated in 1986 Chudnovsky/Chudnovsky:

     K<d,X1,X2,Y1,Y2>:=FieldOfFractions(PolynomialRing(Rationals(),5));
     R<Z1,Z2>:=PolynomialRing(K,2);
     S:=quo<R|X1^3+Y1^3+Z1^3-3*d*X1*Y1*Z1,X2^3+Y2^3+Z2^3-3*d*X2*Y2*Z2>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^3+y1^3+1-3*d*x1*y1);
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^3+y2^3+1-3*d*x2*y2);
     x3:=(y1^2*x2-y2^2*x1)/(x2*y2-x1*y1);
     y3:=(x1^2*y2-x2^2*y1)/(x2*y2-x1*y1);
     S!(x3^3+y3^3+1-3*d*x3*y3);
     // here are the formulas:
     X1Y2:=X1*Y2; X1Z2:=X1*Z2;
     Y1Z2:=Y1*Z2; Y1X2:=Y1*X2;
     Z1X2:=Z1*X2; Z1Y2:=Z1*Y2;
     X3:=Y1X2*Y1Z2-Z1Y2*X1Y2;
     Y3:=X1Z2*X1Y2-Y1X2*Z1X2;
     Z3:=Z1Y2*Z1X2-X1Z2*Y1Z2;
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

2001 Joye Quisquater, 12M + 3add, 7 registers (none preserved):

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

Hessian readdition with X2=1.
2007 Hisil/Carter/Dawson, 5M + 6S + 12add after 3S + 3add + 2times2 to cache Y2^2, Z2^2, 2*Y2*Z2, 2*Y2, 2*Z2:

     K<d,X1,Y2,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),4));
     R<Z1,Z2>:=PolynomialRing(K,2);
     X2:=1;
     S:=quo<R|X1^3+Y1^3+Z1^3-3*d*X1*Y1*Z1,X2^3+Y2^3+Z2^3-3*d*X2*Y2*Z2>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^3+y1^3+1-3*d*x1*y1);
     x2:=X2/Z2; y2:=Y2/Z2;
     S!(x2^3+y2^3+1-3*d*x2*y2);
     x3:=(y1^2*x2-y2^2*x1)/(x2*y2-x1*y1);
     y3:=(x1^2*y2-x2^2*y1)/(x2*y2-x1*y1);
     S!(x3^3+y3^3+1-3*d*x3*y3);
     // cached:
     S0 := Y2^2;
     S1 := Z2^2;
     S2 := (Y2+Z2)^2-S0-S1;
     S3 := 2*Y2;
     S4 := 2*Z2;
     // uncached:
     R0 := X1^2;
     R1 := Y1^2;
     R2 := Z1^2;
     R3 := X1+Y1;
     Y3 := Y1+Z1;
     Y3 := Y3^2;
     Z3 := X1+Z1;
     Z3 := Z3^2;
     X3 := R3^2;
     R3 := Y3-R1;
     R3 := R3-R2;
     Y3 := R0*S2;
     Y3 := Y3-R3;
     R3 := X3-R0;
     R3 := R3-R1;
     Z3 := Z3-R0;
     Z3 := Z3-R2;
     X3 := R1*S4;
     R0 := Z3*S0;
     X3 := X3-R0;
     Z3 := R2*S3;
     R0 := R3*S1;
     Z3 := Z3-R0;
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

Hessian mixed addition.
Sylvester, 10M + 3add, substituting Z2=1 in the above:

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

Hessian addition with Z1=1 and Z2=1.
Sylvester, 8M + 3add:

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

Hessian doubling.
Cauchy-Sylvester doubling law:

     K<d,x1>:=FieldOfFractions(PolynomialRing(Rationals(),2));
     R<y1>:=PolynomialRing(K,1);
     S:=quo<R|x1^3+y1^3+1-3*d*x1*y1>;
     x2:=x1; y2:=y1;
     // the Cauchy-Sylvester doubling law:
     x3:=y1*(1-x1^3)/(x1^3-y1^3);
     y3:=x1*(y1^3-1)/(x1^3-y1^3);
     // map to the Weierstrass curve:
     p1:=12*(d^3-1)/(d+x1+y1); u1:=p1-9*d^2; v1:=3*p1*(y1-x1);
     S!(v1^2-u1^3+27*d*(d^3+8)*u1-54*(d^6-20*d^3-8));
     p2:=12*(d^3-1)/(d+x2+y2); u2:=p2-9*d^2; v2:=3*p2*(y2-x2);
     S!(v2^2-u2^3+27*d*(d^3+8)*u2-54*(d^6-20*d^3-8));
     p3:=12*(d^3-1)/(d+x3+y3); u3:=p3-9*d^2; v3:=3*p3*(y3-x3);
     S!(v3^2-u3^3+27*d*(d^3+8)*u3-54*(d^6-20*d^3-8));
     // add on the Weierstrass curve:
     lambda:=(3*u1^2-27*d*(d^3+8))/(2*v1);
     r3:=lambda^2-u1-u2; s3:=lambda*(u1-r3)-v1;
     // check the answer:
     S!(u3-r3); S!(v3-s3);

Inversion-free version:

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

6M + 3S + 3add, consistent with the "9 multiplications" stated in 1986 Chudnovsky/Chudnovsky:

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

2007 Hisil/Carter/Dawson, 7M + 1S + 8add:

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

2007 Hisil/Carter/Dawson, 7M + 1S + 8add:

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

2007 Hisil/Carter/Dawson, 3M + 6S + 15add + 3times2:

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

2007 Hisil/Carter/Dawson, 3M + 6S + 15add + 3times2:

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

Hessian weakly unified addition.
2001 Joye Quisquater, 12M + 3add, weakly unified, applying the addition formulas to permuted input coordinates:

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

Hessian doubling with Z1=1.
2007 Hisil/Carter/Dawson, with Z1=1, scaled, 3M + 3S + 11add + 3times2:

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

Hessian tripling.
2007 Hisil/Carter/Dawson:

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

2007 Hisil/Carter/Dawson, common subexpressions eliminated, 11M + 4S + 2D + 8add:

     K<d,X1,Y1>:=FieldOfFractions(PolynomialRing(Rationals(),3));
     a:=3*d;
     R<Z1>:=PolynomialRing(K,1);
     X2:=X1; Y2:=Y1; Z2:=Z1;
     S:=quo<R|X1^3+Y1^3+Z1^3-3*d*X1*Y1*Z1>;
     x1:=X1/Z1; y1:=Y1/Z1;
     S!(x1^3+y1^3+1-3*d*x1*y1);
     x2:=y1*(1-x1^3)/(x1^3-y1^3);
     y2:=x1*(y1^3-1)/(x1^3-y1^3);
     S!(x2^3+y2^3+1-3*d*x2*y2);
     x3:=(y1^2*x2-y2^2*x1)/(x2*y2-x1*y1);
     y3:=(x1^2*y2-x2^2*y1)/(x2*y2-x1*y1);
     S!(x3^3+y3^3+1-3*d*x3*y3);
     // here are the formulas:
     XX:=X1^2; A:=XX*X1;
     YY:=Y1^2; B:=YY*Y1;
     ZZ:=Z1^2; C:=ZZ*Z1;
     AB:=A-B;
     BC:=B-C;
     CA:=C-A;
     U:=B*CA;
     V:=A*BC;
     X3:=a*(U*AB-V*BC);
     Y3:=a*(V*AB-U*CA);
     Z3:=(A+B+C)*(BC*CA-AB^2);
     // check:
     S!(x3-X3/Z3); S!(y3-Y3/Z3);

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

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