RR<x,S,t,u0,u1,u2,u3,u4> := PolynomialRing(Integers(),8);
P := S*t^4 + 5; F := x^5 - S*(t^3*x^4 + 2*t^2*x^3 + 2*t*x^2 + x);
I := Ideal([P,F]);
printf("Sanity check that (1+t*x)^5 == 1 mod (P,F).\n"); (1+t*x)^5 -1 in I;
RRbar, pi := quo< RR | I >;
   // this is the C_p group action
D := hom< RR -> RR | x,S,t,u0,
 x*u0 + (1+t*x)*u1,
 x^2*u0 + 2*x*(1+t*x)*u1 + (1+t*x)^2*u2,
 x^3*u0 + 3*x^2*(1+t*x)*u1 + 3*x*(1+t*x)^2*u2 + (1+t*x)^3*u3,
 x^4*u0 + 4*x^3*(1+t*x)*u1 + 6*x^2*(1+t*x)^2*u2 + 4*x*(1+t*x)^3*u3 + (1+t*x)^4*u4 >;
printf("Check that D defines an action of GG_a on PP^4.\n");
[D(D(D(D(D(m))))) - m in Ideal(I) : m in [u0,u1,u2,u3,u4]];

// Start by stating and verifying the answer.
// Prop. Each of the following matrices satisfies
// D(N) = Du*N*(Transpose Du) modulo the relations.
// (The derivation of the Nij is below.)

N34 := -AntisymmetricMatrix([
0,
0,0,
0,0,0,
0,0,0,u0+2*t*u1+t^2*u2
]);

N24 := -AntisymmetricMatrix([
0,
0,0,
0,0,0,
0,0,u0+t*u1,3*u1+3*t*u2
]);

N23 := -AntisymmetricMatrix([
0,
0,0,
0,0,u0,
0,0,4*u1,6*u2
]);

N14 := -AntisymmetricMatrix([
0,
0,0,
0,0,0,
0,u0,2*u1,3*u2
]);

N13 := -AntisymmetricMatrix([
0,
0,0,
0,u0+4*t*u1+6*t^2*u2+4*t^3*u3+t^4*u4,-8*u1-12*t*u2-8*t^2*u3-2*t^3*u4,
0,-16*u1-24*t*u2-16*t^2*u3-4*t^3*u4,-8*S*t^3*u1+48*u2+32*t*u3+8*t^2*u4,+8*S*t^2*u1+8*S*t^3*u2-32*u\
3-8*t*u4
]);

N04 := -AntisymmetricMatrix([
0,
0,0,
0,0,0,
u0+4*t*u1+6*t^2*u2+4*t^3*u3+t^4*u4,-4*u1-6*t*u2-4*t^2*u3-t^3*u4,-S*t^3*u1+6*u2+4*t*u3+t^2*u4,S*t^2\
*u1+S*t^3*u2-4*u3-t*u4
]);

N12 := -AntisymmetricMatrix([
0,
0,u0+3*t*u1+3*t^2*u2+t^3*u3,
0,3*u1+9*t*u2+9*t^2*u3+3*t^3*u4,3*S*t^3*u1-27*u2-21*t*u3-6*t^2*u4,
0,6*S*t^3*u1-54*u2-42*t*u3-12*t^2*u4,-16*S*t^2*u1-24*S*t^3*u2+88*u3+24*t*u4,18*S*t*u1+24*S*t^2*u2+\
18*S*t^3*u3-24*u4
]);

N03 := -AntisymmetricMatrix([
0,
0,0,
u0+3*t*u1+3*t^2*u2+t^3*u3,u1+3*t*u2+3*t^2*u3+t^3*u4,S*t^3*u1-9*u2-7*t*u3-2*t^2*u4,
4*u1+12*t*u2+12*t^2*u3+4*t^3*u4,4*S*t^3*u1-36*u2-28*t*u3-8*t^2*u4,-8*S*t^2*u1-12*S*t^3*u2+44*u3+12\
*t*u4,9*S*t*u1+12*S*t^2*u2+9*S*t^3*u3-12*u4
]);

N02 := -AntisymmetricMatrix([
0,
u0+2*t*u1+t^2*u2,u1+2*t*u2+t^2*u3,
3*u1+6*t*u2+3*t^2*u3,3*u2+6*t*u3+3*t^2*u4,2*S*t^2*u1+4*S*t^3*u2-18*u3-6*t*u4,
6*u2+12*t*u3+6*t^2*u4, 6*S*t^2*u1+12*S*t^3*u2-54*u3-18*t*u4,
 -15*S*t*u1-25*S*t^2*u2-20*S*t^3*u3+30*u4, 18*S*u1+27*S*t*u2+21*S*t^2*u3+6*S*t^3*u4
]);

N01 := -AntisymmetricMatrix([
u0+t*u1,
2*u1+2*t*u2,u2+t*u3,
3*u2+3*t*u3,2*u3+2*t*u4,S*t*u1+2*S*t^2*u2+2*S*t^3*u3-4*u4,
4*u3+4*t*u4,3*S*t*u1+6*S*t^2*u2+6*S*t^3*u3-12*u4,-8*S*u1-14*S*t*u2-12*S*t^2*u3-4*S*t^3*u4,-2*S^2*t\
^3*u1+16*S*u2+13*S*t*u3+ 4*S*t^2*u4
]);

Du := Matrix(5,[1,0,0,0,0, x,(1+t*x),0,0,0, x^2,2*(1+t*x)*x,(1+t*x)^2,0,0, x^3,3*(1+t*x)*x^2,3*(1+t*x)^2*x,(1+t*x)^3,0,
 x^4,4*(1+t*x)*x^3,6*(1+t*x)^2*x^2,4*(1+t*x)^3*x,(1+t*x)^4]);

printf("Check that the 10 matrices Nij have invariance
 Du*N*Transpose(Du) == D(N).\n");
for N in [N01,N02,N03,N04,N12,N13,N14,N23,N24,N34] do
DNN := Matrix(5,[D(m) : m in Eltseq(N)]);
Matrix(5, [pi(m) : m in Eltseq(DNN - Du*N*Transpose(Du))]) eq 0;
end for;

FF5 := FiniteField(5);
Rbar<u0,u1,u2,u3,u4> := PolynomialRing(FF5,5);
bar5 := hom< RR -> Rbar | 0,0,0,u0,u1,u2,u3,u4 >;

Nbar01 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N01)]);
Nbar02 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N02)]);
Nbar03 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N03)]);
Nbar04 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N04)]);
Nbar12 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N12)]);
Nbar13 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N13)]);
Nbar14 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N14)]);
Nbar23 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N23)]);
Nbar24 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N24)]);
Nbar34 := Matrix(Rbar, 5, [bar5(m) : m in Eltseq(N34)]);

printf("The 10 matrices Nij reduced modulo (S,t,p) are:\n");

[Nbar01,Nbar02,Nbar03,Nbar04, Nbar12,Nbar13,Nbar14, Nbar23,Nbar24, Nbar24];

Nbar := &+[1*N : N in
[Nbar01,Nbar02,Nbar03,Nbar04, Nbar12,Nbar13,Nbar14, Nbar23,Nbar24, Nbar34]];

printf("Choose the particular linear combination Nbar01+Nbar03+Nbar24.\n");

Nbar := &+[1*N : N in
[Nbar01, Nbar03, Nbar24]];
Nbar;

printf("Check that Nbar defines a nonsingular elliptic curve over FF5.\n");
E := Scheme(Proj(Rbar), Pfaffians(Nbar,4));
Dimension(E);
IsNonSingular(E);
Genus(Curve(E));