A2<t,u,v>:=AffineSpace(Rationals(),3);

//Surface S1 defined by Equation (2.3)

f:=-16*t^4*u*v + 96*t^2*u*v + t*u^4*v^4 - 6*t*u^4*v^2 + 8*t*u^4*v - 3*t*u^4 -

6*t*u^2*v^4 + 36*t*u^2*v^2 - 48*t*u^2*v + 18*t*u^2 + 8*t*u*v^4 - 48*t*u*v^2

- 64*t*u*v - 24*t*u - 3*t*v^4 + 18*t*v^2 - 24*t*v + 9*t + 48*u*v;

S1:=Scheme(A2,f);

//Surface S2 defined by biquadratic polynomial in u and v in proof of Theorem 2.3

f2:=(16*t^4 - t)*u^2*v^2 + (-16*t^4 - 2*t)*u^2*v + 3*t*u^2 + (-16*t^4 - 2*t)*u*v^2 +

(16*t^4 + 64*t^3 - 4*t - 1)*u*v + (6*t - 3)*u + 3*t*v^2 + (6*t - 3)*v - 9*t

- 9;

S2:=Scheme(A2,f2);

//Construct birational map between S1 and S2

G:=[(u-1)*(v-1)*t+1,u,v];

m1:=map<S1->S2|[(t-1)/((u-1)*(v-1)),u,v]>;

m1inv:=map<S2->S1|G>;

//View f2 as a quadratic in u and use quadratic equation to get another birational surface S3

f3:=(256*t^8 - 128*t^5 + 16*t^2)*v^4 + (-512*t^8 - 2048*t^7 - 128*t^5 - 32*t^4 +

64*t^2 - 8*t)*v^3 + (256*t^8 + 2048*t^7 + 4096*t^6 + 640*t^5 - 64*t^4 -

128*t^3 - 32*t^2 - 40*t + 1)*v^2 + (-384*t^5 + 96*t^4 - 384*t^3 - 192*t^2 -

24*t + 6)*v + 144*t^2 + 72*t + 9-256*u^2;

S3:=Scheme(A2,f3);

//Coefficients of u

a0:=(t^4 - 1/16*t)*v^2 + (-t^4 - 1/8*t)*v + 3/16*t;

b0:=(-t^4 - 1/8*t)*v^2 + (t^4 + 4*t^3 - 1/4*t - 1/16)*v + 3/8*t - 3/16;

c0:=3/16*t*v^2 + (3/8*t - 3/16)*v - 9/16*t - 9/16;

H:=[t,(-b0+u)/(2*a0),v];

m2inv:=map<S3->S2|H>;

_,m2:=IsInvertible(m2inv);

//Use standard equations to get from y^2=quartic (and a point) to y^2=cubic

m3pols:= [

2*t,

(128*t^8*v^2 - 128*t^8*v - 512*t^7*v - 64*t^5*v^2 - 32*t^5*v + 96*t^5 + 128*t^4*u -

8*t^4*v + 168*t^4 + 8*t^2*v^2 + 16*t^2*v - 24*t^2 - 32*t*u - 2*t*v - 6*t)/(t^2 + t

+ 1/4),

(2048*t^12*v^3 - 3072*t^12*v^2 + 1024*t^12*v - 12288*t^11*v^2 + 8192*t^11*v + 16384*t^10*v -

1536*t^9*v^3 + 2304*t^9*v - 768*t^9 + 2048*t^8*u*v - 1024*t^8*u + 2880*t^8*v^2 -

2304*t^8*v + 192*t^8 - 4096*t^7*u - 4608*t^7*v - 768*t^7 + 384*t^6*v^3 + 576*t^6*v^2 -

768*t^6*v - 192*t^6 - 1024*t^5*u*v - 256*t^5*u - 96*t^5*v - 96*t^5 - 64*t^4*u +

132*t^4*v + 204*t^4 - 32*t^3*v^3 - 96*t^3*v^2 + 32*t^3*v + 96*t^3 + 128*t^2*u*v +

128*t^2*u + 12*t^2*v^2 + 40*t^2*v + 12*t^2 - 16*t*u - t*v - 3*t)/(t^3 + 3/2*t^2 +

3/4*t + 1/8)

];

m3invpols:=[

1/2*t,

(1/64*t^21 + 1/4*t^20 - 1/4608*t^19*u - 1/32*t^19 - 1/288*t^18*u - 47/64*t^18 +

1/1327104*t^17*u^2 - 7/512*t^17*u - 41/64*t^17 + 1/82944*t^16*u^2 + 23/3072*t^16*u +

23/32*t^16 + 1/20736*t^15*u^2 + 13/256*t^15*u - 1/576*t^15*v + 57/32*t^15 -

11/663552*t^14*u^2 + 35/1024*t^14*u - 1/576*t^14*v + 3/32*t^14 - 329/663552*t^13*u^2 -

275/4608*t^13*u + 7/9216*t^13*v - 51/32*t^13 - 79/82944*t^12*u^2 - 533/4608*t^12*u +

13/9216*t^12*v - 25/32*t^12 + 1/663552*t^11*u^3 + 43/331776*t^11*u^2 - 1/1327104*t^11*v^2 +

5/3072*t^11*v + 7/16*t^11 + 1/331776*t^10*u^3 + 733/331776*t^10*u^2 + 533/4608*t^10*u -

1/663552*t^10*v^2 + 7/3072*t^10*v + 37/64*t^10 + 1/663552*t^9*u^3 + 3457/1327104*t^9*u^2 +

275/4608*t^9*u - 1/1327104*t^9*v^2 + 1/1024*t^9*v + 1/16*t^9 - 1/110592*t^8*u^3 -

73/165888*t^8*u^2 - 35/1024*t^8*u + 1/221184*t^8*v^2 - 7/3072*t^8*v - 1/8*t^8 -

1/55296*t^7*u^3 - 547/165888*t^7*u^2 - 13/256*t^7*u + 1/110592*t^7*v^2 - 5/3072*t^7*v -

1/32*t^7 - 1/110592*t^6*u^3 - 1327/663552*t^6*u^2 - 23/3072*t^6*u + 1/221184*t^6*v^2 +

5/9216*t^6*v + 1/55296*t^5*u^3 + 7/10368*t^5*u^2 + 7/512*t^5*u - 1/110592*t^5*v^2 +

1/4608*t^5*v + 1/27648*t^4*u^3 + 131/82944*t^4*u^2 + 1/288*t^4*u - 1/55296*t^4*v^2 -

1/4608*t^4*v + 1/55296*t^3*u^3 + 143/331776*t^3*u^2 + 1/4608*t^3*u - 1/110592*t^3*v^2 -

1/4608*t^3*v - 1/82944*t^2*u^3 - 1/2592*t^2*u^2 + 1/165888*t^2*v^2 - 1/41472*t*u^3 -

1/10368*t*u^2 + 1/82944*t*v^2 - 1/82944*u^3 - 1/165888*u^2 + 1/165888*v^2)/(t^17 -

1/72*t^15*u - 2*t^15 - t^14 + 1/20736*t^13*u^2 + 1/72*t^13*u + t^13 + 1/16*t^12*u +

2*t^12 + 1/4*t^11 - 1/2592*t^10*u^2 - 1/16*t^10*u - t^10 - 1/12*t^9*u - 1/2*t^9 +

1/864*t^7*u^2 + 1/12*t^7*u + 1/4*t^7 + 1/36*t^6*u - 1/648*t^4*u^2 - 1/36*t^4*u +

1/1296*t*u^2),

(1/2*t^9 - 1/288*t^7*u - 1/2*t^7 - 1/36*t^6*u + 1/2*t^6 - 1/144*t^4*u - 1/288*t^4*v -

1/2*t^4 - 1/288*t^3*u - 1/288*t^3*v + 1/36*t*u + 1/144*t*v - 1/144*u + 1/144*v)/(t^9 -

1/144*t^7*u - t^7 - 1/2*t^6 + 1/36*t^4*u + 1/2*t^4 - 1/36*t*u)

];

f4:=u^3+(t-1)^2*(((t)^2+8*(t)+1)*u-144*(t)^4)^2-v^2;

S4:=Scheme(A2,f4);

m3:=map<S3->S4|m3pols>;

m3inv:=map<S4->S3|m3invpols>;

m:=m1*m2*m3;

DP:=DefiningPolynomials(m);

K:=FunctionField(A2);

PP<z>:=PolynomialRing(K);

E:=EllipticCurve(z^3+(t-1)^2*(((t)^2+8*(t)+1)*z-144*(t)^4)^2);

//Sections coming from points with u=-3, v=0, and u=0, v=-3

Q1:=E![Evaluate(DP[2],[2*t+1,-3,0]), Evaluate(DP[3],[2*t+1,-3,0])];

Q2:=E![Evaluate(DP[2],[2*t+1,0,-3]), Evaluate(DP[3],[2*t+1,0,-3])];

//3-torsion point

T:=E![0,144*t^4*(t-1)];

//P1 and P2 as in Theorem 2.3

P1:=-Q1-T;

P2:=-Q2+T;

“P_1, P_2, T as in Theorem 2.3”;

P1;P2;

T;

//Construct the solution of (2.3) mentioned before Theorem 2.4

minv:=m3inv*m2inv*m1inv;

DPinv:=DefiningPolynomials(minv);

P3:=T-2*P2;

“Equation for t/u before Theorem 2.4”;

Evaluate(DPinv[1],[t,P3[1],P3[2]])/Evaluate(DPinv[2],[t,P3[1],P3[2]]);

“Equation for u and v before Theorem 2.4”;

Evaluate(DPinv[2],[t,P3[1],P3[2]]);

Evaluate(DPinv[3],[t,P3[1],P3[2]]);