// Construct the schemes A, B and P3B, together with maps between them.
A<s1,s2,s3> := AffineSpace(Rationals(),3);

A4<t1,t2,t3,e> := AffineSpace(Rationals(),4);
B := Scheme(A4, e^4+1);
i := e*e;

phi := map< B->A | [t1^4,t2^4,t3^4] >;

P3<X0,X1,X2,X3> := ProjectiveSpace(Rationals(),3);
P3A4<t1p, t2p, t3p, ep, X0p, X1p, X2p, X3p>, pi := DirectProduct(A4,P3);
P3B := Scheme(P3A4, (e @@ pi[1])^4 + 1);
pi[1] := Restriction(pi[1], P3B, B);

AutBGens := [ iso< B->B | [i*t1,t2,t3,e], [-i*t1,t2,t3,e] >,
              iso< B->B | [t1,i*t2,t3,e], [t1,-i*t2,t3,e] >,
              iso< B->B | [t1,t2,i*t3,e], [t1,t2,-i*t3,e] >,
              iso< B->B | [t1,t2,t3,-e], [t1,t2,t3,-e] >,
              iso< B->B | [t1,t2,t3,e^7], [t1,t2,t3,e^7] > ];

// This reduces a scheme isomorphism to a sensible form.
ReduceIso := func< f | 
  iso< Domain(f) -> Codomain(f) |
    [ CoordinateRing(Domain(f)) ! p : p in DefiningPolynomials(f) ],
    [ CoordinateRing(Codomain(f)) ! p : p in InverseDefiningPolynomials(f) ] 
  >
>;

// Now we can construct the automorphism group of B/A.
// In Magma, things act on the *right*.
AutB, BAction := GenericGroup(AutBGens :
                              Mult := func<f1,f2 | ReduceIso(f1*f2)>);


// Lift an automorphism of B to one of P3B
LiftAut := func< f | 
  map< P3B->P3B |
    [ x @@ pi[1] : x in DefiningPolynomials(f) ] 
      cat [X0p, X1p, X2p, X3p],
    [ x @@ pi[1] : x in InverseDefiningPolynomials(f) ]
      cat [X0p, X1p, X2p, X3p]
  >
>;

// Define the action of AutB on P3B
P3BAction := map< AutB->Parent(IdentityAutomorphism(P3B)) | 
               f :-> LiftAut(BAction(f)) >;

// Construct the 48 lines
L123 := [ Scheme(P3B, [X0p - ep^m*t1p*X1p, t2p*X2p - ep^n*t3p*X3p])
          : n, m in [1,3,5,7] ];
L231 := [ Scheme(P3B, [X0p - ep^m*t2p*X2p, t3p*X3p - ep^n*t1p*X1p])
          : n, m in [1,3,5,7] ];
L312 := [ Scheme(P3B, [X0p - ep^m*t3p*X3p, t1p*X1p - ep^n*t2p*X2p])
          : n, m in [1,3,5,7] ];
L := L123 cat L231 cat L312;

// This is a function for comparing lines.
LineCompare := func<l1,l2 | Dimension(l1 meet l2) eq 4 >;

// Return the range of 16 integers containing i
Range16 := func<i | 
  { ((16 * ((i-1) div 16)) + 1) .. (16 * ((i-1) div 16 + 1)) } >;

// Given an element of AutB, this function returns the permutation
// vector representing its action on the 48 lines
function LAction(f)
  P := [];
  for j := 1 to 48 do
    LL := L[j] @ P3BAction(f);
    _ := exists(P[j]) { i : i in Range16(j) | LineCompare(L[i], LL)};
  end for;
  return P;
end function;

// In fact this is a faithful permutation action, so we'll define
// G to be the associated permutation group.  It's a more convenient
// thing to work with that the FP-group AutB.
G := PermutationGroup<48 | [ LAction(AutB.i) : i in {1..5} ] >;
GAutBIso := hom<G->AutB | [AutB.1, AutB.2, AutB.3, AutB.4, AutB.5] >;
BAction := GAutBIso * BAction;
P3BAction := GAutBIso * P3BAction;

// Create an AutB-module representing this action
Lambda := PermutationModule(G, Integers());

// Return the intersection number of two lines
IntersectionNumber := func<l1,l2 | 
  case< Dimension(l1 meet l2) | 3 : 1, 4 : -2, default : 0 >
>;

// Create a matrix for the intersection form on Lambda
Q := SymmetricMatrix([ IntersectionNumber(L[i],L[j]) : 
                         j in {1 .. i}, i in {1 .. 48} ]);

// Now we want to create the quotient G-module of Lambda by the kernel
// of Q.  But we can't do that.  First we create the quotient of the
// underlying Z-module.
LambdaZ := VectorSpace(Lambda); // poorly-named function
PicZ, PicZMap := quo< LambdaZ | Kernel(Q) >;

// To compute the G-action on the quotient, we need a matrix representing
// the quotient map PicZMap...
PicZMatrix := Matrix( [ LambdaZ.i @ PicZMap : i in {1..48} ] );

// ... and a section (in this case a *right* inverse) of it.
PicZSection := Solution(PicZMatrix, IdentityMatrix(Integers(), 20));

// Now we can form the quotient G-module
Pic := GModule(G, [ PicZSection * M * PicZMatrix : 
                       M in ActionGenerators(Lambda) ] );

// This function takes rationals [a0,a1,a2,a3] and finds the associated
// subgroup of G.  We take the point a in A and lift it to b in B; then
// split b into prime components and take the stabiliser of one of them
// in G.
FindGroupSlow := function(a)
  a := Cluster(A ! [a[2]/a[1], a[3]/a[1], a[4]/a[1]]);
  b := a @@ phi;
  bc := { Cluster(B, DefiningPolynomials(x)) : x in PrimaryComponents(b) };
  bcset := GSet(G, bc, map< car<bc,G>->bc | x :-> x[1] @ BAction(x[2]) > );
  return Stabiliser(G, bcset, Rep(bcset));
end function;

// Solve sum n_i e_i = y in an Abelian group A
AbSolve := function(A, e, y)
  // Get a matrix of relations for A, that is, such that A is the 
  // quotient of Z^n by the lattice generated by the rows of A.
  M := Matrix( [ Eltseq(LHS(r) - RHS(r)) : r in Relations(A) ] );
  // Get a matrix representing the elements e_i
  N := Matrix( [ Eltseq(x) : x in e ] );
  // Combine these matrices, and solve the resulting linear system
  v := Solution( VerticalJoin(N, M), Vector(Eltseq(y)) );
  // Throw away the bits corresponding to the relations, and 
  // return the rest.
  return ColumnSubmatrix(v, NumberOfRows(N))[1];
end function;

SimplifyMap := func<m | 
  hom< Domain(m) -> Codomain(m) | 
    [ m(Domain(m).i) : i in [1 .. Ngens(Domain(m))] ] > >;

DirectSumPower := function(G,n)
  H := G;
  i := [ IdentityHomomorphism(G) ];
  p := [ IdentityHomomorphism(G) ];
  for k := 2 to n do
    H, j1, j2, q1, q2 := DirectSum(H, G);
    i := [ SimplifyMap(f * j1) : f in i ] cat [ j2 ];
    p := [ SimplifyMap(q1 * f) : f in p ] cat [ q2 ];
  end for;
  return H, i, p;
end function;

Hom2 := function(G,H)
  X, t := Hom(G,H);
  n := Ngens(G);
  m := Ngens(X);
  Hn,i,p := DirectSumPower(H,n);
  u := function(f) 
    a := AbSolve(Hn,	
      [ &+ [ i[j](G.j @ t(X.k)) : j in [1..n] ] : k in [1..m] ],
      &+ [ i[j](G.j @ f) : j in [1..n] ] );
    return &+ [ a[j] * X.j : j in [1..m] ];
  end function;
  return X, t, u;
end function;

HomFunctor := function(f,C)
  AC, tA, u := Hom2(Domain(f), C);
  BC, tB := Hom(Codomain(f), C);
  return AC, BC, hom<BC->AC | [ u(f * tB(BC.i)) : i in [1..Ngens(BC)] ] >,
    tA, tB;
end function;

// This does the same as FindGroupSlow, but in a faster and less 
// natural way.  It also works over other fields.
GetRoot := func< n, S, r | [ ( (p eq -1) select (1 - Sign(n))/2
                                         else Valuation(n, p) ) 
                    * (4 div r) : p in S ] >;

FindGroup := function(V, f)
  // Find the primes we're interested in
  S := [-1] cat V`DQBadPrimes;

  X := AbelianGroup([4,4,4,2,4]);
  Q4 := V`DQFourthRoots;

  a := V`DQCoefficients;
  g := hom<X->Q4 | 
         [ Q4 ! GetRoot(x,S,4) : x in [ a[i]/a[1] : i in {2..4} ]] 
         cat [ Q4 ! GetRoot(2,S,2), Q4 ! GetRoot(-4,S,4) ] > * f;

  C4 := CyclicGroup(GrpAb, 4);

  Y := Codomain(g);
  XDual, tX, uX := Hom2(X, C4);
  YDual, tY, uY := Hom2(Y, C4);
  XXDual, tXX, uXX := Hom2(sub<X | X.5>, C4);

  gdual := hom< YDual->XDual | 
    [ uX(g * tY(YDual.i)) : i in [1..Ngens(YDual)] ] >; 
  idual := hom< XDual->XXDual |
    [ uXX(tX(XDual.i)) : i in [1..Ngens(XDual)] ]>;

  HH := Image(gdual);

  v := [1,1,1,2,1];
  p := map<XDual -> G | 
    h :-> &* [ G.i ^ (Eltseq(X.i @ tX(h))[1] div v[i]) : i in [5..1 by -1] ] >;

  Z := { h : h in Generators(HH) | X.5 @ tX(h) eq 3*C4.1 } join
        { -h : h in Generators(HH) | X.5 @ tX(h) eq C4.1 };

  g5 := IsEmpty(Z) select Identity(G) else p(Rep(Z));

  return sub<G | { p(h) : h in Generators(HH meet Kernel(idual)) }, g5 >;
end function;

LocalFourthRoots := function(V, p)
  n := (p eq 2) select 16 else p;
  U, f := UnitGroup(Integers(n));
  case (p mod 4):
    when 1: // Already have a fourth root of unity
      U1 := U;
      m := map<Integers(n) -> U1 | x :-> x @@ f>;
    when 3: // Adjoin a fourth root of unity
      U1 := CyclicGroup(GrpAb, 2*(p-1));
      m := map<Integers(n) -> U1 | x :-> U1 ! Eltseq(x @@ f)>;
    when 2: // Special case
      U1 := AbelianGroup([4,4]);
      m := map<Integers(n) -> U1 | x :-> U1 ! Eltseq(x @@ f)>;
  end case;
  pZ := FreeAbelianGroup(1);
  A, i1, i2 := DirectSum(pZ, U1);
  A4, q := quo< A | [ 4*g : g in Generators(A) ] >;
  return hom< V`DQFourthRoots -> A4 | 
    [ (x eq p) select q(i1(pZ.1)) else q(i2(m(x))) : 
      x in [-1] cat V`DQBadPrimes ] >;
end function;  

AddAttribute(Sch, "DQBadPrimes");
AddAttribute(Sch, "DQCoefficients");
AddAttribute(Sch, "DQGlobalGaloisGroup");
AddAttribute(Sch, "DQFourthRoots");
AddAttribute(Sch, "DQLocalGaloisGroups");
AddAttribute(Sch, "DQCohomologyModule");
AddAttribute(Sch, "DQH1Pic");

DiagonalQuarticSurface := function(a)
  V := Scheme(P3, a[1]*X0^4 + a[2]*X1^4 + a[3]*X2^4 + a[4]*X3^4);
  V`DQCoefficients := a;

  // Get bad primes
  V`DQBadPrimes := Setseq( &join { Seqset(PrimeDivisors(a[i])) : i in {1..4} } 
                            join { 2 } );
  Sort(~V`DQBadPrimes);

  // Get the group generated by fourth roots of the bad primes, and e
  V`DQFourthRoots := AbelianGroup([4] cat [ 4 : p in V`DQBadPrimes ]);

  // Initialise a few other things
  V`DQLocalGaloisGroups := {};

  return V;
end function;

BadPrimes := func<V | V`DQBadPrimes>;

GaloisGroup := function(V, p)
  if not assigned V`DQGlobalGaloisGroup then
    V`DQGlobalGaloisGroup := FindGroup(V, 
      IdentityHomomorphism(V`DQFourthRoots));
  end if;
  if p eq 0 then
    return V`DQGlobalGaloisGroup;
  else
    if not exists(H){ t[2] : t in V`DQLocalGaloisGroups | t[1] eq p } then
      H := FindGroup(V, LocalFourthRoots(V,p));
      Include(~V`DQLocalGaloisGroups, <p,H>);
    end if;
    return H;
  end if;
end function;

PicCM := function(V)
  H := GaloisGroup(V, 0);
  if not assigned V`DQCohomologyModule then
    V`DQCohomologyModule := CohomologyModule(H, Restriction(Pic, H));
  end if;
  return V`DQCohomologyModule;
end function;

H1Pic := function(V)
  if not assigned V`DQH1Pic then
    V`DQH1Pic := CohomologyGroup(PicCM(V), 1);
  end if;
  return V`DQH1Pic;
end function;

RestrictionKernel := function(V, p)
  M := PicCM(V);
  X := H1Pic(V);
  c := [ OneCocycle(M, X.i) : i in [ 1 .. #Moduli(X) ] ];
  H := GaloisGroup(V,p);
  if Order(H) ne 1 then
    N := CohomologyModule(H, Restriction(Pic, H));
    Y := CohomologyGroup(N,1);
    res := hom< X -> Y | [ IdentifyOneCocycle(N, f) : f in c ] >;
    return Kernel(res);
  end if;
  return X;
end function;