printf("==== 1. Derivation of ring of invariants for TO5 over B ====\n");
printf("==== 2. Reduction mod (S,t,5) to Godeaux quintic hypersurfaces
over FF5 and proof of nonsingularity ====\n");

// Derivation of ring of invariants for TO5 acting on
// k[u0,u1,u2,u3,u4]. The invariants are as for mu5
// acting by 1/5(0,1,2,3,4), mostly familiar from Godeaux
// quintics.

// Reduction mod (S,t,5) to quintics of PP^4. Proof that the
// general image of this reduction defines a nonsingular
// quotient Calabi-Yau 3-fold Y5/al5 and Godeaux surface
// T5/al5 with al5 contained in Pic^0.

RR<S,t,x,u0,u1,u2,u3,u4> := PolynomialRing(Integers(),8);
P := 5 + S*t^4;
F := x^5 - S*(t^3*x^4 + 2*t^2*x^3 + 2*t*x^2 + x);
T := 1+t*x;
IPF := Ideal([P,F]);
T^5-1 in IPF;
RRbar, pi := quo< RR | [P,F] >;

f := hom< RR -> RR | S, t, x, u0, x*u0 + T*u1,
 x^2*u0 + 2*x*T*u1 + T^2*u2,
 x^3*u0 + 3*x^2*T*u1 + 3*x*T^2*u2 + T^3*u3,
 x^4*u0 + 4*x^3*T*u1 + 6*x^2*T^2*u2 + 4*x*T^3*u3 + T^4*u4 >;

function Mult(T, q)
  // multiplicity of prime factor q in term T. This assumes lots
  // of things, but I only use it for T = integer * monomial
XX := Factorization(T);
   if &and[q ne m[1] : m in XX] then return 0;
   else return [m : m in XX | m[1] eq q][1,2];
   end if;
end function;

function vt5(F)  // pseudo-valuation at t and 5, with v(5) = 4.
  return Min([Mult(T, t) + 4*Mult(T, 5) : T in Terms(F)]);
end function;

function psi(X) X51 := [T : T in Terms(X) | Mult(T,5) ge 1];
while X51 ne [] do X := X - (&+X51 div 5)*P;
X51 := [T : T in Terms(X) | Mult(T,5) ge 1];
end while;
return X;
end function;

v0 := u0;
v1 := u0 + t*u1;
v2 := u0 + 2*t*u1 + t^2*u2;
v3 := u0 + 3*t*u1 + 3*t^2*u2 + t^3*u3;
v4 := u0 + 4*t*u1 + 6*t^2*u2 + 4*t^3*u3 + t^4*u4;

printf("==== Invariants have two homogeneity properties, with =====\n");
printf("== deg ui = 1, deg S,t = 0, wt ui = i, wt t = -1, wt S = 4: ===\n");
printf("== deg vi = 1, wt vi = 0, so wt = power of t cancelled ===\n");
printf("====\n");

printf("==== u0 is invariant, deg 1, wt 0 =====\n");
X := u0;
f(X) eq X;

printf("==== B1 from v1*v4 / t^2 gives u2, deg 2, wt 2 =====\n");
X := v1*v4 - u0^2;
B1 := psi(X) div t^2;  B1;
f(B1) - B1 in Ideal([P,F]);

printf("==== C1 from v1^2*v3 / t^3, gives u3, deg 3, wt 3 =====\n");
X := v1^2*v3 - 3*u0*v1*v4 + 2*u0*u0^2;
C1 := psi(X) div t^3; C1;
f(C1) - C1 in Ideal([P,F]);

printf("==== B2 from v2*v3 / t^4 gives u4, deg 2, wt 4 =====\n");
X := v1*v4 + v2*v3 - 2*u0^2;
B2 := psi(X) div t^4; B2;
f(B2) - B2 in Ideal([P,F]);

printf("==== C2 from v1*v2^2 / t^5 gives S*u1 or u1*u2^2, deg 3, wt 5 =====\n");
X := v1*v2^2 + 2*v1^2*v3 + 8*u0*v2*v3 - 11*u0^3;
C2 := psi(X) div t^5; C2;
f(C2) - C2 in Ideal([P,F]);

printf("==== C3 from v2*v4^2 / t^6, deg 3, wt 6 =====\n");
X := v2*v4^2 + v1^2*v3 - 16*u0*v2*v3 - 17*u0*v1*v4 + 31*u0^3;
C3 := psi(X) div t^6; C3;
f(C3) - C3 in Ideal([P,F]);

printf("==== D1 from v1^3*v2 / t^7, deg 4, wt 7 =====\n");
X := v1^3*v2 + 6*v1*v2*v3*v4 + 39*v1^2*v4^2 - 4*u0*v1^2*v3
 - u0*v1*v2^2 - 3*u0*v2*v4^2 - 3*u0^2*v1*v4 + 23*u0^2*v2*v3 - 58*u0^4;
D1 := psi(X) div t^7; D1;
f(D1) - D1 in Ideal([P,F]);

printf("==== C4 from v3^2*v4 / t^8, deg 3, wt 8 =====\n");
X := v3^2*v4 + 6*v2*v4^2 - 14*v1*v2^2 - 24*v1^2*v3 - 128*u0*v2*v3
 - 13*u0*v1*v4 + 172*u0^3;
C4 := psi(X) div t^8; C4;
f(C4) - C4 in Ideal([P,F]);

printf("==== D2 from v1*v3^3 / t^9, deg 4, wt 9 =====\n");
X := 679*u0^4 - 88*u0^2*v1*v4 - 431*u0^2*v2*v3 - 32*u0*v1^2*v3 - 25*u0*v1*v2^2 +
    27*u0*v2*v4^2 + 2*u0*v3^2*v4 - 2*v1^3*v2 - 112*v1^2*v4^2 - 19*v1*v2*v3*v4 +
    v1*v3^3;
D2 := psi(X) div t^9; D2;
f(D2) - D2 in Ideal([P,F]);

printf("==== D3 from v2^3*v4 / t^10, deg 4, wt 10 =====\n");
X := v1*v3^3 - 2*v1*v2*v3*v4 - 872*v1^2*v4^2 + v2^3*v4 + 6*v2^2*v3^2 -
33*u0*v1^2*v3 - 24*u0*v1*v2^2 + 17*u0*v2*v4^2 + u0*v3^2*v4 +
1585*u0^2*v1*v4 - 57*u0^2*v2*v3 - 623*u0^4;
D3 := psi(X) div t^10; D3;
f(D3) - D3 in Ideal([P,F]);

printf("==== D4 from v3*v4^3 / t^12, deg 4, wt 12 =====\n");
X := v1*v3^3 + v2^3*v4 + 24*v2^2*v3^2 + v3*v4^3 + v1^3*v2 +
4499*v1^2*v4^2 + 26*v1*v2*v3*v4 + 48*u0*v1^2*v3 + 23*u0*v1*v2^2 -
2*u0*v2*v4^2 - 2*u0*v3^2*v4 - 9147*u0^2*v1*v4 + 253*u0^2*v2*v3 + 4274*u0^4;
D4 := psi(X) div t^12; D4;
f(D4) - D4 in Ideal([P,F]);

printf("==== E1 from v1^5 / t^5, deg 5, wt 5 =====\n");
X := v1^5 - v0^5;
E1 := psi(X) div t^5; E1;
f(E1) - E1 in Ideal([P,F]);

printf("==== E2 from v2^5 / t^10, deg 5, wt 10 =====\n");
X := v2^5 - 10*u0^2*v1*v2^2 - 10*u0*v1^3*v2 - 7*v1^5 + 26*u0^5;
E2 := psi(X) div t^10; E2;
f(E2) - E2 in Ideal([P,F]);

printf("==== E3 from v3^5 / t^15, deg 5, wt 15 =====\n");
X := v3^5 - 60*u0*v1*v3^3 - 270*u0*v2^2*v3^2 - 810*v1^2*v2*v3^2
     - 1620*v1*v2^3*v3 - 243*v2^5 + 15060*u0^3*v2*v3
     + 67770*u0^2*v1^2*v3 + 203310*u0^2*v1*v2^2 + 406620*u0*v1^3*v2
     + 60993*v1^5 - 750751*u0^5;
E3 := psi(X) div t^15; E3;
f(E3) - E3 in Ideal([P,F]);

printf("==== E4 from v4^5 / t^20, deg 5, wt 20 =====\n");
X := v4^5 - 80*u0*v3*v4^3 - 480*v1*v2*v4^3 - 1920*v1*v3^2*v4^2
 - 4320*v2^2*v3*v4^2 - 7680*v2*v3^3*v4 - 1024*v3^5
 + 45180*u0^2*v2*v4^2 + 120480*u0^2*v3^2*v4 + 120480*u0*v1^2*v4^2
 + 2891520*u0*v1*v2*v3*v4 + 1285120*u0*v1*v3^3 + 1084320*u0*v2^3*v4
 + 4337280*u0*v2^2*v3^2 + 1285120*v1^3*v3*v4 + 4337280*v1^2*v2^2*v4
 + 11566080*v1^2*v2*v3^2 + 17349120*v1*v2^3*v3 + 1951776*v2^5
 - 60060080*u0^3*v1*v4 - 360360480*u0^3*v2*v3 - 1441441920*u0^2*v1^2*v3
 - 3243244320*u0^2*v1*v2^2 - 5765767680*u0*v1^3*v2 - 768769024*v1^5
 + 11593285251*u0^5;
E4 := psi(X) div t^20; E4;
f(E4) - E4 in Ideal([P,F]);

Inv5 := [u0^5, u0^3*B1, u0^3*B2, u0^2*C1, u0^2*C2, u0^2*C3, u0^2*C4,
 u0*B1^2, u0*B1*B2, u0*B2^2, u0*D1, u0*D2, u0*D3, u0*D4,
 B1*C1, B1*C2, B1*C3, B1*C4, B2*C1, B2*C2, B2*C3, B2*C4,
 E1,E2,E3,E4]; #Inv5;

// reduction modulo (S,t,p,x,u0) to case of Godeaux surface.
// The Calabi-Yau 3-fold case is similar, but the calculation
// is bigger.

printf("\n==== Godeaux quintics F5 in PP^3 over FF5 ====\n");
printf("==== choice [ 17, 22 ] of r(Inv5) ====\n");

FF5 := FiniteField(5);
RRbar<u1,u2,u3,u4> := PolynomialRing(FF5,4);
r := hom< RR -> RRbar | 0,0,0,0,u1,u2,u3,u4 >;
[r(m) : m in Inv5][15..26];

f := u1^5+u2^5+u3^5+u4^5 + 1*r(Inv5[17]) + 1*r(Inv5[22]);
Y := Scheme(Proj(RRbar), f);
printf("==== Dimension, degree and reduced degree ====\n");
Dimension(SingularSubscheme(Y));
Degree(SingularSubscheme(Y));
Degree(ReducedSubscheme(SingularSubscheme(Y)));
f;

printf("==== Choice [ 21, 22 ], differs by ====\n");
r(Inv5[21] - Inv5[17]);

f := u1^5+u2^5+u3^5+u4^5 + 1*r(Inv5[21]) + 1*r(Inv5[22]);
Y := Scheme(Proj(RRbar), f);
Dimension(SingularSubscheme(Y));
Degree(SingularSubscheme(Y));
Degree(ReducedSubscheme(SingularSubscheme(Y)));
f;

printf("\n==== CY quintics Y5 in PP^4 over FF5 ====\n");
printf("==== Sum of invariants [1,14,17,23,24,25,26] of r(Inv5) ====\n");
printf("==== Restrict mod (S,t,x) to FF5, then work in affine piece ui = 1 ====\n");
printf("==== Check SingularScheme is a union of 40 disjoint orbits of al5 ====\n");

// Passing to affine pieces is a little kludge to deal with the fact
// that Magma's ReducedSubscheme does not work properly in PP^4.

RRbar<u0,u1,u2,u3,u4> := PolynomialRing(FF5,5);
r := hom< RR -> RRbar | 0,0,0,u0,u1,u2,u3,u4 >;
L := [
    u0^5,
    u0^4*u2 + 4*u0^3*u1^2,
    u0^4*u4 + u0^3*u1*u3 + 3*u0^3*u2^2,
    4*u0^4*u3 + 3*u0^3*u1*u2 + 3*u0^2*u1^3,
    3*u0^3*u2*u3 + 2*u0^2*u1^2*u3 + u0^2*u1*u2^2,
    4*u0^3*u2*u4 + u0^3*u3^2 + u0^2*u1^2*u4 + 3*u0^2*u1*u2*u3 + u0^2*u2^3,
    u0^3*u4^2 + 2*u0^2*u1*u3*u4 + u0^2*u2^2*u4 + u0^2*u2*u3^2,
    u0^3*u2^2 + 3*u0^2*u1^2*u2 + u0*u1^4,
    u0^3*u2*u4 + 4*u0^2*u1^2*u4 + u0^2*u1*u2*u3 + 3*u0^2*u2^3 + 4*u0*u1^3*u3 +
        2*u0*u1^2*u2^2,
    u0^3*u4^2 + 2*u0^2*u1*u3*u4 + u0^2*u2^2*u4 + u0*u1^2*u3^2 + u0*u1*u2^2*u3 +
        4*u0*u2^4,
    4*u0^3*u3*u4 + 3*u0^2*u1*u2*u4 + 4*u0^2*u1*u3^2 + 4*u0^2*u2^2*u3 +
        3*u0*u1^3*u4 + 4*u0*u1^2*u2*u3 + 3*u0*u1*u2^3,
    4*u0^2*u2*u3*u4 + 4*u0^2*u3^3 + u0*u1^2*u3*u4 + 3*u0*u1*u2^2*u4 +
        3*u0*u1*u2*u3^2,
    2*u0^2*u2*u4^2 + u0^2*u3^2*u4 + 3*u0*u1^2*u4^2 + 3*u0*u1*u2*u3*u4 +
        u0*u1*u3^3 + u0*u2^3*u4 + u0*u2^2*u3^2,
    u0^2*u4^3 + 3*u0*u1*u3*u4^2 + 4*u0*u2^2*u4^2 + 3*u0*u2*u3^2*u4 + 4*u0*u3^4,
    4*u0^3*u2*u3 + u0^2*u1^2*u3 + 3*u0^2*u1*u2^2 + 2*u1^5,
    3*u0^2*u2^2*u3 + 4*u0*u1^2*u2*u3 + u0*u1*u2^3 + 3*u1^4*u3 + 4*u1^3*u2^2,
    4*u0^2*u2^2*u4 + u0^2*u2*u3^2 + 2*u0*u1^2*u2*u4 + 4*u0*u1^2*u3^2 +
        3*u0*u1*u2^2*u3 + u0*u2^4 + 4*u1^4*u4 + 2*u1^3*u2*u3 + 4*u1^2*u2^3,
    u0^2*u2*u4^2 + 4*u0*u1^2*u4^2 + 2*u0*u1*u2*u3*u4 + u0*u2^3*u4 + u0*u2^2*u3^2
        + 3*u1^3*u3*u4 + 4*u1^2*u2^2*u4 + 4*u1^2*u2*u3^2,
    4*u0^3*u3*u4 + 3*u0^2*u1*u2*u4 + 4*u0^2*u1*u3^2 + 2*u0^2*u2^2*u3 +
        3*u0*u1^3*u4 + 3*u0*u1^2*u2*u3 + 4*u0*u1*u2^3 + 3*u1^4*u3 + 4*u1^3*u2^2,
    3*u0^2*u2*u3*u4 + 2*u0*u1^2*u3*u4 + u0*u1*u2^2*u4 + 3*u0*u1*u2*u3^2 +
        4*u0*u2^3*u3 + 2*u1^3*u3^2 + 2*u1^2*u2^2*u3 + 3*u1*u2^4,
    4*u0^2*u2*u4^2 + u0^2*u3^2*u4 + u0*u1^2*u4^2 + 2*u0*u1*u2*u3*u4 + u0*u1*u3^3
        + 3*u0*u2^3*u4 + 3*u0*u2^2*u3^2 + u1^3*u3*u4 + 3*u1^2*u2^2*u4 +
        3*u1^2*u2*u3^2 + 3*u2^5,
    u0^2*u4^3 + 3*u0*u1*u3*u4^2 + 4*u0*u2^2*u4^2 + u0*u2*u3^2*u4 +
        2*u1^2*u3^2*u4 + 2*u1*u2^2*u3*u4 + u1*u2*u3^3 + 3*u2^4*u4 + 3*u2^3*u3^2,
    u1^5,
    u2^5,
    u3^5,
    u4^5
];
L eq [r(m) : m in Inv5];

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

f := &+[s(Inv5[i]) : i in [1,14,17,23,24,25,26]]; // Should give nonsingular X
Y := Scheme(Spec(RRa), f);
Dimension(SingularSubscheme(Y));
Degree(SingularSubscheme(Y));
Degree(ReducedSubscheme(SingularSubscheme(Y)));
f;