RR<x,S,t,u0,u1,u2,u3,u4,u5,u6> := PolynomialRing(Integers(),10);
u := [u0,u1,u2,u3,u4,u5,u6]; // basis of regular representation Vreg.
P := S*t^6 + 7; F := -x^7 + S*(t^5*x^6+3*t^4*x^5+5*t^3*x^4+5*t^2*x^3+3*t*x^2 + x);
IPF := Ideal([P,F]);
T := 1+t*x; T^7-1 in IPF;
RRbar,pi := quo< RR | IPF >;
 printf("Now over B1 = ZZ[S,t]/(S*t^6+7). Check D defines an action of GG_7 on PP^7.\n");
Du := Matrix(7,[1,0,0,0,0,0,0,
 x,T,0,0,0,0,0,
 x^2,2*T*x,T^2,0,0,0,0,
 x^3,3*T*x^2,3*T^2*x,T^3,0,0,0,
 x^4,4*T*x^3,6*T^2*x^2,4*T^3*x,T^4,0,0,
 x^5,5*T*x^4,10*T^2*x^3,10*T^3*x^2,5*T^4*x,T^5,0,
 x^6,6*T*x^5,15*T^2*x^4,20*T^3*x^3,15*T^4*x^2,6*T^5*x,T^6]);
// Du;
&and[m in IPF : m in Eltseq(Du^7-1)];
Du*Matrix(1,u);
D := hom< RR -> RR | [x,S,t] cat Eltseq(Du*Matrix(1,u)) >;
[D(D(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;
v5 := u0 + 5*t*u1 + 10*t^2*u2 + 10*t^3*u3 + 5*t^4*u4 + t^5*u5; D(v5) eq T^5*v5;
v6 := u0 + 6*t*u1 + 15*t^2*u2 + 20*t^3*u3 + 15*t^4*u4 + 6*t^5*u5 + t^6*u6; D(v6) eq T^6*v6;
// Construct B : v1 -> v3, so u1 -> Imu1 = 3*u1 + 3*t*u2 + t^2*u3
Imu1 := (v3 - u0) div t; // Imu1;
Imu2 := (v6 - u0 - 2*t*Imu1) div t^2;
Imu3 := (v2 - u0 - 3*t*Imu1 - 3*t^2*Imu2 + P*(5*t^2*u2 + t*u1)) div t^3;
Imu4 := (v5 - u0 - 4*t*Imu1 - 6*t^2*Imu2 - 4*t^3*Imu3 + P*(t*u1 + 8*t^2*u2 - 18*t^3*u3)) div t^4;
Imu5 := (v1 - v0 - 5*t*Imu1 - 10*t^2*Imu2 - 10*t^3*Imu3 - 5*t^4*Imu4
     + P*(2*t*u1 + 15*t^2*u2 - 55*t^3*u3 + 25*t^4*u4)) div t^5;
Imu6 := (v4 - v0 - 6*t*Imu1 - 15*t^2*Imu2 - 20*t^3*Imu3 - 15*t^4*Imu4 - 6*t^5*Imu5
     + P*(2*t*u1 + 21*t^2*u2 - 124*t^3*u3 + 107*t^4*u4 - 15*t^5*u5)) div t^6;
printf("======= Images B(ui) =========\n");
[Imu1, Imu2, Imu3, Imu4, Imu5, Imu6];
B := hom< RR -> RR | x,S,t,u0,Imu1, Imu2, Imu3, Imu4, Imu5, Imu6 >;
&and[B(B(B(B(B(B(RR.i)))))) - RR.i in IPF : i in [1..10]];
Bu := Matrix([[&+([0] cat [x div n : x in Terms(B(m)) | IsDivisibleBy(x,n)]) : n in u] : m in u]);
// Bu;
printf("======= Check Bu^6 = 1 and Bu*Du = D^3*Bu =========\n");
&and[m in IPF : m in Eltseq(Bu^6 - 1)];
&and[m in IPF : m in Eltseq(Bu*Du-Du^3*Bu)];
// [[i,j] : j in [i+1..6], i in [0..5]]; // remember how to use Magma 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..6]];
b := [&+([0] cat [T div u[i+1] : T in Terms(w2) | IsDivisibleBy(T, u[i+1])]) : i in [0..6]];
return
&+([0] cat [a[i+1]*b[j+1]*u[i+1]*u[j+1] : j in [i+1..6], i in [0..5]])
- &+([0] cat [a[i+1]*b[j+1]*u[i+1]*u[j+1] : j in [0..i-1], i in [1..6]]);
end function;
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..6], k in [0..5]]
  : j in [i+1..6], i in [0..5]]);
printf("======= Matrix Dw too big to print. Check D^7 = 1. =========\n");
// Dw;
Matrix(21, [pi(m) : m in Eltseq(Dw^7)]) eq 1;
Bw := Matrix([
[&+ ([0] cat [T div (u[k+1]*u[l+1]) : T in Terms(wedge(B(u[i+1]),B(u[j+1])))
  | IsDivisibleBy(T, u[k+1]*u[l+1]) ])
  : l in [k+1..6], k in [0..5]]
  : j in [i+1..6], i in [0..5]]);
// Bw;
printf("======= Check matrix Bw has Bw^6 = 1 and Bw*Dw = Dw^3*Bw. =========\n");
&and[m in IPF : m in Eltseq(Bw^6 - 1)];
&and[m in IPF : m in Eltseq(Bw*Dw-Dw^3*Bw)];

RRM<m010,m011,m012,m013,m014,m015,m016,
m020,m021,m022,m023,m024,m025,m026,
m030,m031,m032,m033,m034,m035,m036,
m040,m041,m042,m043,m044,m045,m046,
m050,m051,m052,m053,m054,m055,m056,
m060,m061,m062,m063,m064,m065,m066,
m120,m121,m122,m123,m124,m125,m126,
m130,m131,m132,m133,m134,m135,m136,
m140,m141,m142,m143,m144,m145,m146,
m150,m151,m152,m153,m154,m155,m156,
m160,m161,m162,m163,m164,m165,m166,
m230,m231,m232,m233,m234,m235,m236,
m240,m241,m242,m243,m244,m245,m246,
m250,m251,m252,m253,m254,m255,m256,
m260,m261,m262,m263,m264,m265,m266,
m340,m341,m342,m343,m344,m345,m346,
m350,m351,m352,m353,m354,m355,m356,
m360,m361,m362,m363,m364,m365,m366,
m450,m451,m452,m453,m454,m455,m456,
m460,m461,m462,m463,m464,m465,m466,
m560,m561,m562,m563,m564,m565,m566> := PolynomialRing(RR,147);
[RRM.i : i in [29..35]];
[RRM.i : i in [36..42]];
[RRM.i : i in [43..49]];
[RRM.i : i in [50..56]];
// RRM;
incl := hom< RR -> RRM | x,S,t,u0,u1,u2,u3,u4,u5,u6 >;
DuM := Matrix(7, [incl(x) : x in Eltseq(Du)]);
DwM := Matrix(21, [incl(x) : x in Eltseq(Dw)]);
MM := Matrix(7,[RRM.i : i in [1..147]]);
X := Eltseq(DwM*MM - MM*DuM);   // the condition for MM to be Cp-invariant
#X; // X[1..5];

L := [[i,j] : j in [i+1..7], i in [1..6]]; // lists the 21 [i,j] in order

function MS(SS); // SS is a linear form in 147 mijk.
  // MS translates it to the skew matrix with ijth entry = sum SSijk u^k
MSS := ZeroMatrix(RR,7,7);
for I in [1..21] do MSS[L[I][1], L[I][2]] :=
    &+([0] cat [MonomialCoefficient(SS,RRM.(7*(I-1)+k))*u[k] : k in [1..7]]); end for;
return MSS - Transpose(MSS);
end function;

NN := Matrix([[RR!pi(MonomialCoefficient(xx, RRM.i)) : i in [1..147]] : xx in X]);
printf("==========  Solution with m56 has u0 ===========\n");
Si := [141..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
S56 := Sol[1,1];
MS56 := MS(S56); // MS56

printf("==========  Solution with m46 has u0 ===========\n");
Si := [134..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
S46 := Sol[1,1];
MS46 := MS(S46); // MS46

printf("==========  Solution with m45 has u0 ===========\n");
Si := [127..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
S45 := Sol[1,1];
MS45 := MS(S45); // MS45;

printf("==========  Solution with m36 has u0 ===========\n");
Si := [120..126] cat [134..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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;
S36 := Sol[1,1];
MS36 := MS(S36); // MS36;

printf("==========  Solution with m35 has u0 ===========\n");
Si := [113..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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;
S35 := Sol[1,1];
MS35 := MS(S35); // MS35;

printf("==========  Solution with m26 has u0 ===========\n");
Si := [99..105] cat [120..126] cat [134..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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;
S26 := Sol[1,1];
MS26 := MS(S26); // MS26;

printf("==========  Solution with m34 has u0 ===========\n");
Si := [106..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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;
S34 := Sol[1,1];
MS34 := MS(S34); // MS34;

printf("==========  Solution with m25 has u0 ===========\n");
Si := [92..105] cat [113..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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;
S25 := Sol[1,1];
MS25 := MS(S25); // MS25;

printf("==========  Solution with m16 has u0 ===========\n");
Si := [71..77] cat [99..105] cat [120..126] cat [134..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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;
S16 := Sol[1,1];
MS16 := MS(S16); // MS16;

printf("==== To proceed, we need new matrix and new function ===========\n");
NN := Matrix([[MonomialCoefficient(xx, RRM.i) : i in [1..147]] : xx in X]);
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("==========  Solution with m24 has u0 ===========\n");
Si := [85..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
Y := Sol[1,1];
Not(Y);
Y1 := Y + 75*S56; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 75*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 45*S36 + 20*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 15*S26 - 15*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 5*S25 + 3*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S24 := Y; S24;
MS24 := MS(S24); // MS24;

printf("==========  Solution with m15 has u0 ===========\n");
Si := [64..77] cat [92..105] cat [113..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
Y := Sol[1,1];
Not(Y);
Y1 := Y + 24*S56; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 24*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 18*S36 + 4*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 12*S26 - 3*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 6*S16 + 2*S25; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S15 := Y;
MS15 := MS(S15); // MS15;

printf("==========  Solution with m06 has u0 ===========\n");
Si := [36..42] cat [71..77] cat [99..105] cat [120..126] cat [134..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
Y := Sol[1,1]; 
Not(Y);
Y1 := Y+S56; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y-S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y+S36; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y-S26; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y+S16; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S06 := Y;
MS06 := MS(S06); // MS06;

printf("==========  Solution with m23 has u0 ===========\n");
Si := [78..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[2,1]; 
Not(Y);
Y1 := Y - 50*S56; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 60*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 45*S36 - 20*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 20*S26 + 20*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 10*S25 - 6*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 4*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S23 := Y;
MS23 := MS(S23); // MS23;

printf("==========  Solution with m14 has u0 ===========\n");
Si := [57..77] cat [85..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[2,1];
Not(Y);
Y1 := Y + 9*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 15*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 15*S26 + 15*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 15*S16 - 5*S25 - 6*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 3*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S14 := Y; 
MS14 := MS(S14); // MS14;

printf("==========  Solution with m05 has u0 ===========\n");
Si := [29..42] cat [64..77] cat [92..105] cat [113..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[2,1]; Not(Y);
Y1 + S46;
Y1 := Y + S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - S36 - S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + S26 + S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - S16 - S25; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + S06 + S15; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S05 := Y;
MS05 := MS(S05); // MS05;

printf("==========  Solution with m13 has u0 ===========\n");
Si := [50..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[5,1]; Not(Y);
Y1 := Y - 8*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 9*S36 + 10*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 15*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 10*S16 + 10*S25 + 8*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 8*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 2*S14 + 3*S23; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S13 := Y;
MS13 := MS(S13); // MS13;

printf("==========  Solution with m04 has u0 ===========\n");
Si := [22..42] cat [57..77] cat [85..105] cat [106..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[6,1]; Not(Y);
Y1 := Y - 2*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 15*S36 + 14*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 45*S26 - 15*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 75*S16 + 5*S25 + 7*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 105*S06 + 5*S15 - 5*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 15*S05 + 3*S14; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S04 := Y;
MS04 := MS(S04); // MS04;

printf("==========  Solution with m12 has u0 ===========\n");
Si := [43..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[6,1]; Not(Y); // -306*m360 - 265*m450
Y1 := Y + 306*S36 + 265*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 816*S26 - 675*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 1020*S16 + 1005*S25 + 399*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 925*S15 - 640*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 535*S14 + 219*S23; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 174*S13; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S12 := Y;
MS12 := MS(S12); // MS12;

printf("==========  Solution with m03 has u0 ===========\n");
Si := [8..42] cat [50..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[10,1]; Not(Y);
Y1 := Y - 4*S36 + 2*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 34*S26 + 10*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 339*S16 - 61*S25 + 18*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 1444*S06 + 56*S15 + 4*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 349*S05 - 61*S14 + 17*S23; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 42*S04 + 12*S13; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S03 := Y;
MS03 := MS(S03); // MS03;

printf("==========  Solution with m02 has u0 ===========\n");
Si := [8..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[15,1]; Not(Y);
Y1 := Y + 5*S56; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 9*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 12*S36 + 4*S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 14*S26 + 7*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 15*S16 + 9*S25 + 3*S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 15*S06 - 10*S15 - 5*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 10*S05 + 6*S14 + 2*S23; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 6*S04 - 3*S13; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 3*S03 + S12; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S02 := Y + 2*S*S35;
MS02 := MS(S02); // MS02;

printf("==========  Solution with m01 has u0 ===========\n");
Si := [1..147]; RRM.(Si[1]);
NN4 := Submatrix(NN,[1..147], 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]);
[Terms(Sol[n,1])[1] : n in [1..NumberOfRows(Sol)]];
Y := Sol[16,1]; Not(Y);
Y1 := Y - S56; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 2*S46; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 3*S36 - S45; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 4*S26 + 2*S35; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 5*S16 - 3*S25 - S34; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 6*S06 + 4*S15 + 2*S24; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 5*S05 - 3*S14 - S23; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 4*S04 + 2*S13; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y - 3*S03 - S12; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
Y1 := Y + 2*S02; Y := (Y1 - P*Not(Y1) div 7) div t; Not(Y);
S01 := Y;
MS01 := MS(S01); // MS01;

printf("\n==== Check that the 21 matrices Nij have invariance
 Du*N*Transpose(Du) == D(N). ====\n");
for N in [MS01,MS02,MS03,MS04,MS05,MS06, MS12,MS13,MS14,MS15,MS16,
MS23,MS24,MS25,MS26, MS34,MS35,MS36, MS45,MS46,MS56] do
DNN := Matrix(7,[D(m) : m in Eltseq(N)]);
Matrix(7, [pi(m) : m in Eltseq(DNN - Du*N*Transpose(Du))]) eq 0;
end for;

printf("\n===== Nonsingularity checks. =====\n");
cf := [ 5, 4, 4, 5, 3, 6, 5, 5, 1, 6, 5, 6, 2, 1, 3, 4, 4, 4, 6, 3, 3 ];
SS := &+[cf[i]*SSS[i] : i in [1..21]] where SSS is [S01, S01, S02, S03, S04, S05, S06, 
S12, S13, S14, S15, S16, 
S23, S24, S25, S26, 
S34, S35, S36, 
S45, S46, S56]; // SS;

N := MS(SS);

printf("=== Check that the CY Y^3_14 in PP^6 / FF7 has 14 x A6 points. ===\n");
FF7 := FiniteField(7);
RRbar<u0,u1,u2,u3,u4,u5,u6> := PolynomialRing(FF7,7);
r := hom< RR -> RRbar | 0,0,0,u0,u1,u2,u3,u4,u5,u6 >;

IPf := Ideal([r(m) : m in Pfaffians(N,6)]);

Y := Scheme(Proj(RRbar), IPf);
Dimension(SingularSubscheme(Y));
Degree(SingularSubscheme(Y));

RRa<u0,u1,u3,u4,u5,u6> := PolynomialRing(FF7,6);  // Restrict to affine piece u2 = 1.
u2 := 1;
s := hom< RR -> RRa | 0,0,0,u0,u1,u2,u3,u4,u5,u6>;

IPfa := Ideal([s(m) : m in Pfaffians(N,6)]);
Ya := Scheme(Spec(RRa), IPfa);
Dimension(SingularSubscheme(Ya));
Degree(SingularSubscheme(Ya));
Degree(ReducedSubscheme(SingularSubscheme(Ya)));

printf("Check that the surface T^3_14 in PP^5 / FF7 has 10 A6 points.\n");
FF7 := FiniteField(7);
RRbar<u1,u2,u3,u4,u5,u6> := PolynomialRing(FF7,6);
r := hom< RR -> RRbar | 0,0,0,0,u1,u2,u3,u4,u5,u6 >;

IPf := Ideal([r(m) : m in Pfaffians(N,6)]);

Y := Scheme(Proj(RRbar), IPf);
Dimension(SingularSubscheme(Y));
Degree(SingularSubscheme(Y));

RRa<u1,u2,u3,u4,u6> := PolynomialRing(FF7,5);  // Restrict to affine piece u2 = 1.
u5 := 1;
s := hom< RR -> RRa | 0,0,0,0,u1,u2,u3,u4,u5,u6>;

IPfa := Ideal([s(m) : m in Pfaffians(N,6)]);
Ya := Scheme(Spec(RRa), IPfa);
Dimension(SingularSubscheme(Ya));
Degree(SingularSubscheme(Ya));
Degree(ReducedSubscheme(SingularSubscheme(Ya)));