// The derivation

RR<x,S,t,u0,u1,u2,u3,u4> := PolynomialRing(Integers(),8);
u := [u0,u1,u2,u3,u4]; // basis of regular representation Vreg.
P := S*t^4 + 5; F := -x^5 + S*(t^3*x^4+2*t^2*x^3+2*t*x^2 + x);
IPF := Ideal([P,F]);
T := 1+t*x; T^5-1 in IPF; // sanity check to ensure P, F are right
RRbar,pi := quo< RR | IPF >;
 printf("Now over B1 = ZZ[S,t]/(S*t^4+5). Check D defines an action of
TO_5 on PP^4.\n");
// The matrix Du is the group action on Vreg.
Du := Matrix(5,[1,0,0,0,0,
 x,T,0,0,0,
 x^2,2*T*x,T^2,0,0,
 x^3,3*T*x^2,3*T^2*x,T^3,0,
 x^4,4*T*x^3,6*T^2*x^2,4*T^3*x,T^4]);
Du;
&and[m in IPF : m in Eltseq(Du^5-1)];
D := hom< RR -> RR | [x,S,t] cat Eltseq(Du*Matrix(1,u)) >;
[D(D(D(D(D(m))))) - m in IPF : m in u];
printf("======= Check the vi are (1+t*x)^i eigenforms =========\n");
v0 := u0;
v1 := u0 + t*u1; D(v1) eq T*v1;
v2 := u0 + 2*t*u1 + t^2*u2; D(v2) eq T^2*v2;
v3 := u0 + 3*t*u1 + 3*t^2*u2 + t^3*u3; D(v3) eq T^3*v3;
v4 := u0 + 4*t*u1 + 6*t^2*u2 + 4*t^3*u3 + t^4*u4; D(v4) eq T^4*v4;
// [[i,j] : j in [i+1..6], i in [0..5]]; // Remember how Magma orders lists

function wedge(w1,w2)
a := [&+([0] cat [T div u[i+1] : T in Terms(w1) | IsDivisibleBy(T,
u[i+1])]) : i in [0..4]];
b := [&+([0] cat [T div u[i+1] : T in Terms(w2) | IsDivisibleBy(T,
u[i+1])]) : i in [0..4]];
return
&+([0] cat [a[i+1]*b[j+1]*u[i+1]*u[j+1] : j in [i+1..4], i in [0..3]])
- &+([0] cat [a[i+1]*b[j+1]*u[i+1]*u[j+1] : j in [0..i-1], i in [1..4]]);
end function;
// The matrix Dw gives group action on wedge^2 Vreg.
Dw := Matrix([
[&+ ([0] cat [T div (u[k+1]*u[l+1]) : T in Terms(wedge(D(u[i+1]),D(u[j+1])))
  | IsDivisibleBy(T, u[k+1]*u[l+1]) ])
  : l in [k+1..4], k in [0..3]]
  : j in [i+1..4], i in [0..3]]);
printf("======= Matrix Dw. Check Dw^5 = 1. =========\n");
// Dw;
Matrix(10, [pi(m) : m in Eltseq(Dw^5)]) eq 1;

// Introduce 10 x 5 new variables for Hom(Vreg, wedge2 Vreg
RRM<m010,m011,m012,m013,m014,
m020,m021,m022,m023,m024,
m030,m031,m032,m033,m034,
m040,m041,m042,m043,m044,
m120,m121,m122,m123,m124,
m130,m131,m132,m133,m134,
m140,m141,m142,m143,m144,
m230,m231,m232,m233,m234,
m240,m241,m242,m243,m244,
m340,m341,m342,m343,m344> := PolynomialRing(RR,50);
RRM;

incl := hom< RR -> RRM | x,S,t,u0,u1,u2,u3,u4 >;
DuM := Matrix(5, [incl(x) : x in Eltseq(Du)]);
DwM := Matrix(10, [incl(x) : x in Eltseq(Dw)]);
MM := Matrix(5,[RRM.i : i in [1..50]]);
X := Eltseq(DwM*MM - MM*DuM);   // the condition for MM to be TOp-invariant
printf("There are 50 linear conditions for MM to be TOp-invariant. The first 5 are:\n");

#X;
X[1..5];

NN := Matrix([[RR!pi(MonomialCoefficient(xx, RRM.i)) : i in [1..50]] :
xx in X]);
printf("======= Solve B*M = M*A for M with m34 has u0  ===========\n");
Si := [46..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol[1,1];
S34 := Sol[1,1];

printf("======= Solve B*M = M*A for M with m24 has u0  ===========\n");
Si := [41..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol[1,1];
S24 := Sol[1,1];

printf("======= Solve B*M = M*A for M with m23 has u0  ===========\n");
Si := [36..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol;
S23 := Sol[1,1];

printf("======= Solve B*M = M*A for M with m14 has u0  ===========\n");
Si := [31..35] cat [41..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol;
S14 := Sol[1,1];
[S14,S23,S24,S34];

NN := Matrix([[(MonomialCoefficient(xx, RRM.i)) : i in [1..50]] : xx in X]);
printf("======= New matrix, new function No_t in background ===========\n");

// Sum of terms in Y not divisible by t.
function Not(Y); temp := 0;
for T in Terms(Y) do F := Factorisation(T);
  if &and[ not F[i,1] eq t : i in [1..#F-1]] then temp := temp + T;
end if; end for; return temp;
end function;

printf("======= Solve B*M = M*A for M with m13 has u0  ===========\n");
Si := [26..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
// Sol;
Y := Sol[1,1];
Not(Y);  // Take out m560 by hand by subtracting a multiple of S56
Y1 := Y + 8*S34;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y - 8*S24;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y + 4*S14 + 2*S23;
Y := (Y1 - P*Not(Y1) div 5) div t;
Y;
S13 := Y;

printf("======= Solve B*M = M*A for M with m04 has u0  ===========\n");
Si := [16..20] cat [31..35] cat [41..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
// Sol;
Y := Sol[1,1];
Not(Y);
Y1 := Y+S34;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y-S24;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y+S14;
Y := (Y1 - P*Not(Y1) div 5) div t;
Y;
S04 := Y;

printf("======= Solve B*M = M*A for M with m12 has u0  ===========\n");
Si := [21..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol;
Y := Sol[3,1];
Not(Y);
Y1 := Y-6*S34;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y+8*S24;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y-6*S14-3*S23;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y+3*S13;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y;
S12 := Y;

printf("======= Solve B*M = M*A for M with m03 has u0  ===========\n");
Si := [11..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol;
Y := Sol[6,1];
Not(Y);
Y1 := Y-3*S34;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y+4*S24;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y-4*S14-S23;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y+4*S04+S13;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y;
S03 := Y;

printf("======= Solve B*M = M*A for M with m02 has u0  ===========\n");
Si := [6..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol;
Y := Sol[3,1]; Y;
Not(Y);
Y1 := Y + 3*S14 - S23;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y - 9*S04;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
Y1 := Y + 3*S03;
Y := (Y1 - P*Not(Y1) div 5) div t;
Not(Y);
S02 := Y;
Y;

printf("======= Solve B*M = M*A for M with m01 has u0  ===========\n");
Si := [1..50]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..50], Si); NumberOfColumns(NN4); Rank(NN4);
Sol := Matrix(RRM, #Si, [incl(m) : m in
Eltseq(Matrix(MinimalBasis(Kernel(Transpose(NN4)))))]) *
Matrix(1,[RRM.i : i in Si]);
Sol;
Y := Sol[7,1];
Y;
Not(Y);
Y1 := Y - S34; Y := (Y1 - P*Not(Y1) div 5) div t; Not(Y);
Y1 := Y + 2*S24; Y := (Y1 - P*Not(Y1) div 5) div t; Not(Y);
Y1 := Y - 3*S14 - S23; Y := (Y1 - P*Not(Y1) div 5) div t; Not(Y);
Y1 := Y + 4*S04 + 2*S13; Y := (Y1 - P*Not(Y1) div 5) div t; Not(Y);
Y1 := Y - 3*S03 - S12; Y := (Y1 - P*Not(Y1) div 5) div t; Not(Y);
Y1 := Y + 2*S02; Y := (Y1 - P*Not(Y1) div 5) div t; Not(Y);
S01 := Y;
Y;

// Given SS, for example SS := S12, do
Lij := [ [1,2], [1,3], [1,4], [1,5], [2,3], [2,4], [2,5], [3,4], [3,5], [4,5]];

SS := S12;
M12 := ZeroMatrix(RR,5,5);
for I in [1..10] do M12[Lij[I][1], Lij[I][2]] :=
    &+([0] cat [MonomialCoefficient(SS,RRM.(5*(I-1)+k))*u[k] : k in [1..5]]); end for;
M12 := M12 - Transpose(M12);

S12;
M12;

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 >;

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]);

for N in [M12] do
DNN := Matrix(5,[D(m) : m in Eltseq(N)]);
Matrix(5, [pi(m) : m in Eltseq(DNN)])
eq Matrix(5, [pi(m) : m in Eltseq(Du*N*Transpose(Du))]);
end for;

for SS in [S01,S02,S03,S04, S12,S13,S14, S23,S24, S34] do
  Mat := ZeroMatrix(RR,5,5);
for I in [1..10] do Mat[Lij[I][1], Lij[I][2]] :=
    &+([0] cat [MonomialCoefficient(SS,RRM.(5*(I-1)+k))*u[k] : k in [1..5]]); end for;
Mat := Mat - Transpose(Mat); Mat;
DMat := Matrix(5,[D(m) : m in Eltseq(Mat)]);
Matrix(5, [pi(m) : m in Eltseq(DMat)])
eq Matrix(5, [pi(m) : m in Eltseq(Du*Mat*Transpose(Du))]);
end for;