// Computations for "Lines on Fermat surfaces"
// Very inefficiently implemented. For section 7, see checkall.mg

Z:=Integers();

m:=4;
b2:=22;
rho:=20;
q:=3;
Fg<om>:=ext<GF(3)|Polynomial([2,1,1])>;
g:=om^2;

P3<x,y,z,w>:=ProjectiveSpace(Fg,3);
S:=Scheme(P3,x^m+y^m+z^m+w^m);

function GlobalLine(j,k,l)
  if j eq 1 then return Scheme(P3,[om*g^k*x-y,om*g^l*z-w]); end if;
  if j eq 2 then return Scheme(P3,[om*g^k*x-z,om*g^l*y-w]); end if;
  if j eq 3 then return Scheme(P3,[om*g^k*x-w,om*g^l*y-z]); end if;
end function;

GlobalLines:=[GlobalLine(j,k,l): j in [1,2,3], k in [0..m-2], l in [0..m-3]] cat
             [GlobalLine(1,m-1,0)] cat [GlobalLine(2,0,2)];

M:=Matrix(rho,rho,[1+Dimension(l1 meet l2) : l1,l2 in GlobalLines]);
for i in [1..rho] do M[i,i]:=2-m; end for;

assert Determinant(M) eq -64;


lp:=Scheme(P3,[y-x-w,z-x+w]);
lpp:=Scheme(P3,[y/g-x-w,z-x+w]);


alldivs:=GlobalLines cat [lp,lpp];
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in alldivs]);
for i in [1..rho] do N[i,i]:=2-m; end for;
for j in [rho+1..b2] do N[j,j]:=2-m; end for;
Rank(N);
det:=Determinant(N);
if det ne 0 then Factorization(det); end if;



// =====================================




m:=5;
b2:=(m-1)*(m^2-3*m+3)+1;
rho:=3*(m-1)*(m-2)+1;
l:=b2-rho;
q:=4;
Fg<g>:=ext<GF(2)|Polynomial([1,1,1,1,1])>;

P3<x,y,z,w>:=ProjectiveSpace(Fg,3);
S:=Scheme(P3,x^m+y^m+z^m+w^m);

function GlobalLine(j,k,l)
  if j eq 1 then return Scheme(P3,[g^k*x+y,g^l*z+w]); end if;
  if j eq 2 then return Scheme(P3,[g^k*x+z,g^l*y+w]); end if;
  if j eq 3 then return Scheme(P3,[g^k*x+w,g^l*y+z]); end if;
end function;

GlobalLines:=[GlobalLine(j,k,l): j in [1,2,3], k in [0..m-2], l in [1..m-2]] cat
             [GlobalLine(1,m-1,1)];

M:=Matrix(rho,rho,[1+Dimension(l1 meet l2) : l1,l2 in GlobalLines]);
for i in [1..rho] do M[i,i]:=2-m; end for;

assert Determinant(M) eq 5^12;

alpha:=g^3+g^2+1;
lp:=Scheme(P3,[y-alpha*x-alpha^2*w,z-alpha^2*x+alpha*w]);
Dp:=lp;

function Action(a)
  a0,a1,a2:=Explode(a);
  f:=map<P3->P3 | [g^a0*x,g^a1*y,g^a2*z,w]>;
  return f;
end function;

table:=[32..39] cat [44] cat [80..84] cat [93,95];

B2:=[Action([j,k,l])(Dp) 
       where j is Z!((nu-5*k-l-1)/25)
       where k is Z!(((nu-l-1) mod m^2)/5)
       where l is (nu-1) mod m : nu in table];

alldivs:=GlobalLines cat B2;
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in alldivs]);
for i in [1..rho] do N[i,i]:=2-m; end for;
for j in [rho+1..b2] do N[j,j]:=2-m; end for;
Rank(N);
det:=Determinant(N);
if det ne 0 then Factorization(det); end if;



// =====================================


m:=7;
b2:=(m-1)*(m^2-3*m+3)+1;
rho:=3*(m-1)*(m-2)+1;
l:=b2-rho;
p:=13;
F:=GF(p);
Fg<g>:=ext<F|Polynomial([1,5,1])>;

P3<x,y,z,w>:=ProjectiveSpace(Fg,3);
S:=Scheme(P3,x^m+y^m+z^m+w^m);
S2:=Scheme(P3,x^(2*m)+y^(2*m)+z^(2*m)+w^(2*m));
phi:=map<P3->P3 | [x^2,y^2,z^2,w^2]>;

function GlobalLine(j,k,l)
  if j eq 1 then return Scheme(P3,[g^k*x+y,g^l*z+w]); end if;
  if j eq 2 then return Scheme(P3,[g^k*x+z,g^l*y+w]); end if;
  if j eq 3 then return Scheme(P3,[g^k*x+w,g^l*y+z]); end if;
end function;

GlobalLines:=[GlobalLine(j,k,l): j in [1,2,3], k in [0..m-2], l in [1..m-2]] cat
             [GlobalLine(1,m-1,1)];

M:=Matrix(rho,rho,[1+Dimension(l1 meet l2) : l1,l2 in GlobalLines]);
for i in [1..rho] do M[i,i]:=2-m; end for;

assert Determinant(M) eq 7^48;

alpha:=2;
beta:=3*(g-4);
lp:=Scheme(P3,[y-alpha*x-beta*w,z-alpha*w-beta*x]);
Dp:=phi(lp);

function Action(a)
  a0,a1,a2:=Explode(a);
  f:=map<P3->P3 | [g^a0*x,g^a1*y,g^a2*z,w]>;
  return f;
end function;

table:=[[j,k,l] : j in [0..m-2], k in [0..m-2], l in [1..m-2]] cat 
       [[j,0,0] : j in [0..m-2]] cat [[m-1,m-2,m-2]];

function nu(a)
  j,k,l:=Explode(a);
  if k eq 0 and l eq 0 then return b2-m+1+j; end if; 
  if j eq m-1 and k eq m-2 and l eq m-2 then return b2; end if; 
  return 1+j+(m-1)*k+(m-1)^2*(l-1);
end function;

for r in [1..b2] do 
  assert nu(table[r]) eq r;
end for;

Bp:=[Action(table[nu])(Dp) : nu in [1..b2]];
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in Bp]);
for j in [1..b2] do N[j,j]:=-8; end for;
Rank(N);
det:=Determinant(N);
if det ne 0 then Factorization(det); end if;
assert det eq 2^38*7^2*13^48;

Bpprime:=[Action(table[(31*nu mod b2)])(Dp) : nu in [1..l]];
alldivs:=GlobalLines cat Bpprime;
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in alldivs]);
for i in [1..rho] do N[i,i]:=2-m; end for;
for j in [rho+1..b2] do N[j,j]:=-8; end for;
Rank(N);
det:=Determinant(N);
if det ne 0 then Factorization(det); end if;


// =====================================

m:=11;
b2:=(m-1)*(m^2-3*m+3)+1;
rho:=3*(m-1)*(m-2)+1;
l:=b2-rho;
q:=2^5;
f:=CyclotomicPolynomial(11);
Fg<g>:=ext<GF(2)|f>;

P3<x,y,z,w>:=ProjectiveSpace(Fg,3);
S:=Scheme(P3,x^m+y^m+z^m+w^m);
S3:=Scheme(P3,x^(3*m)+y^(3*m)+z^(3*m)+w^(3*m));
phi:=map<P3->P3 | [x^3,y^3,z^3,w^3]>;

function GlobalLine(j,k,l)
  if j eq 1 then return Scheme(P3,[g^k*x+y,g^l*z+w]); end if;
  if j eq 2 then return Scheme(P3,[g^k*x+z,g^l*y+w]); end if;
  if j eq 3 then return Scheme(P3,[g^k*x+w,g^l*y+z]); end if;
end function;

GlobalLines:=[GlobalLine(j,k,l): j in [1,2,3], k in [0..m-2], l in [1..m-2]] cat
             [GlobalLine(1,m-1,1)];

M:=Matrix(rho,rho,[1+Dimension(l1 meet l2) : l1,l2 in GlobalLines]);
for i in [1..rho] do M[i,i]:=2-m; end for;
assert Determinant(M) eq 11^192;

alpha:=g^8+g^7+g^6+g^5+g^4+g^3;
beta:=alpha+1;
lp:=Scheme(P3,[y-alpha*x-beta*w,z-alpha*w-beta*x]);
Dp:=phi(lp);
assert Dp subset S;

table:=[[j,k,l] : j in [0..m-2], k in [0..m-2], l in [1..m-2]] cat 
       [[j,0,0] : j in [0..m-2]] cat [[m-1,m-2,m-2]];

function Action(a)
  a0,a1,a2:=Explode(a);
  f:=map<P3->P3 | [x, g^a0*y,g^a1*z,g^a2*w]>;
  return f;
end function;

function nu(a)
  j,k,l:=Explode(a);
  if k eq 0 and l eq 0 then return b2-m+1+j; end if; 
  if j eq m-1 and k eq m-2 and l eq m-2 then return b2; end if; 
  return 1+j+(m-1)*k+(m-1)^2*(l-1);
end function;

for r in [1..b2] do 
  assert nu(table[r]) eq r;
end for;

Bp:=[Action(table[nu])(Dp) : nu in [1..b2]];
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in Bp]);
for j in [1..b2] do N[j,j]:=-23; end for;
Rank(N);
det:=Determinant(N);
if det ne 0 then Factorization(det); end if;

/*
output:
[ <2, 1200>, <3, 2>, <11, 2>, <23, 64>, <43, 24>, <67, 8>, <131, 16>, <197, 4>, <307,
8>, <331, 8>, <463, 12>, <593, 8>, <3541, 8> ]
*/

Bpprime:=[Action(table[(253*nu mod b2)])(Dp) : nu in [1..l]];
alldivs:=GlobalLines cat Bpprime;
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in alldivs]);
for i in [1..rho] do N[i,i]:=2-m; end for;
for j in [rho+1..b2] do N[j,j]:=-23; end for;
Rank(N);
det:=Determinant(N);
if det ne 0 then Factorization(det); end if;

/*
output:
[ <2, 1202>, <5, 4>, <7, 4>, <23, 48>, <43, 16>, <131, 16>, <439, 2> ]
*/

/*
output for 311 instead of 253:
[ <2, 1202>, <7, 4>, <19, 2>, <23, 48>, <43, 16>, <73, 2>, <131, 16>, <198109, 2>, <1809553, 2> ]
*/


// =====================================

m:=13;
b2:=(m-1)*(m^2-3*m+3)+1;
rho:=3*(m-1)*(m-2)+1;
l:=b2-rho;
q:=5^2;
Fg<g>:=ext<GF(5)|Polynomial([1,1,-1,1,1])>;

P3<x,y,z,w>:=ProjectiveSpace(Fg,3);
S:=Scheme(P3,x^m+y^m+z^m+w^m);
S3:=Scheme(P3,x^(2*m)+y^(2*m)+z^(2*m)+w^(2*m));
phi:=map<P3->P3 | [x^2,y^2,z^2,w^2]>;

function GlobalLine(j,k,l)
  if j eq 1 then return Scheme(P3,[g^k*x+y,g^l*z+w]); end if;
  if j eq 2 then return Scheme(P3,[g^k*x+z,g^l*y+w]); end if;
  if j eq 3 then return Scheme(P3,[g^k*x+w,g^l*y+z]); end if;
end function;

GlobalLines:=[GlobalLine(j,k,l): j in [1,2,3], k in [0..m-2], l in [1..m-2]] cat
             [GlobalLine(1,m-1,1)];

M:=Matrix(rho,rho,[1+Dimension(l1 meet l2) : l1,l2 in GlobalLines]);
for i in [1..rho] do M[i,i]:=2-m; end for;
assert Determinant(M) eq 13^300;

alpha:=2*g^3 + 2*g^2 + g;
beta:=-g^2 + -g + 3;
lp:=Scheme(P3,[y-alpha*x-beta*w,z-alpha*w-beta*x]);
Dp:=phi(lp);
assert Dp subset S;

table:=[[j,k,l] : j in [0..m-2], k in [0..m-2], l in [1..m-2]] cat 
       [[j,0,0] : j in [0..m-2]] cat [[m-1,m-2,m-2]];

function Action(a)
  a0,a1,a2:=Explode(a);
  f:=map<P3->P3 | [x, g^a0*y,g^a1*z,g^a2*w]>;
  return f;
end function;

function nu(a)
  j,k,l:=Explode(a);
  if k eq 0 and l eq 0 then return b2-m+1+j; end if; 
  if j eq m-1 and k eq m-2 and l eq m-2 then return b2; end if; 
  return 1+j+(m-1)*k+(m-1)^2*(l-1);
end function;

for r in [1..b2] do 
  assert nu(table[r]) eq r;
end for;

Bp:=[Action(table[nu])(Dp) : nu in [1..b2]];
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in Bp]);
for j in [1..b2] do N[j,j]:=-20; end for;
Rank(N);
det:=Determinant(N);

if det ne 0 then Factorization(det); end if;

/*
output:
[ <2, 26>, <3, 192>, <5, 912>, <13, 2>, <53, 24>, <79, 24>, <103, 32>, <181, 8>, <233, 8>, <313, 8>, <677, 16>, <883, 4>, <2003, 8>, <2729, 8>, <3847, 8> ]
*/

Bpprime:=[Action(table[(5*nu mod b2)])(Dp) : nu in [1..l]];
alldivs:=GlobalLines cat Bpprime;
N:=Matrix(b2,b2,[(1+Dimension(l1 meet l2))*(Degree(l1 meet l2)) : l1,l2 in alldivs]);
for i in [1..rho] do N[i,i]:=2-m; end for;
for j in [rho+1..b2] do N[j,j]:=-20; end for;
Rank(N);
det:=Determinant(N);

if det ne 0 then Factorization(det); end if;

/*
output:
[ <2, 4>, <3, 144>, <5, 912>, <53, 16>, <103, 32>, <677, 16>, <1151, 2>, <40627, 2>, <42702482453593, 2>, <247634616308749, 2> ]
*/