// This file contains MAGMA code implementing the algorithm
// described in "Computing Galois groups by specialisation".
//
// The main function in FindGroup, at the bottom of the file.
// LocalFourthRoots is also important.  The other functions are
// utilities to implement some functionality not already present
// in Magma.


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


// Turn a map which is stored as a composition of several components
// into a single component
SimplifyMap := func<m | 
  hom< Domain(m) -> Codomain(m) | 
    [ m(Domain(m).i) : i in [1 .. Ngens(Domain(m))] ] > >;

// Form G + ... + G (n times)
// return the sum together with n injections and n projections.
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;


// Return Hom(G,H) as an abstract group together with two functions:  
// t takes an element of the group and returns a map from G to H;
// u takes a map from G to H and returns an element of the group.
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;


// Return the map
//    \sqrt[4]{Q_S} / Q_S  -->  \sqrt[4]{Q_p} / Q_p*
// for the prime p.
// If p=0, return the identity map.
LocalFourthRoots := function(p, Q4, S)
  if p eq 0 then
    return IdentityHomomorphism(Q4);
  end if;
  // Choose n such that ( Z_p* / (Z_p*)^4 ) = ( Z/nZ* / (Z/nZ*)^4 )
  n := (p eq 2) select 16 else p;
  // Get Z/nZ* and the reduction map (well, its inverse actually)
  U, f := UnitGroup(Integers(n));
  // Extend by adjoining a fourth root of unity, if necessary, to
  // get \sqrt[4]{Z/nZ*}
  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;
  // Take the direct sum with the free group generated by p^(1/4)
  pZ := FreeAbelianGroup(1);
  A, i1, i2 := DirectSum(pZ, U1);
  // Compute A / 4A
  A4, q := quo< A | [ 4*g : g in Generators(A) ] >;
  // Compute the map from Q4 to A4
  return hom< Q4 -> A4 | 
    [ (x eq p) select q(i1(pZ.1)) else q(i2(m(x))) : x in S ] >;
end function;  


// Return the element of \sqrt[4]{Q_S} / Q_S corresponding to the
// rth root of n.  (n in Q_S, r = 2 or 4)
GetRoot := func< n, S, r | [ ( (p eq -1) select (1 - Sign(n))/2
                                         else Valuation(n, p) ) 
                    * (4 div r) : p in S ] >;

// Define the group G.
// Magma wants us to put the generators in this order, so we have to.
G := PolycyclicGroup<t,s4,s3,s2,s1 |
  t^2, s4^2, s3^4, s2^4, s1^4,
  s4^t = s4, s3^t = s3^3, s2^t = s2^3, s1^t = s1^3,
  s3^s4 = s3, s2^s4 = s2, s1^s4 = s1,
  s2^s3 = s2, s1^s3 = s1, s1^s2 = s1 >;

// This function finds the Galois group of the extension of Q generated
// by an eighth root of unity and the fourth roots of all the a_i
// (i = 1..3).
//
// Arguments:
//   a = sequence of three rational numbers
//   p = prime (for local group) or 0 (for global group). 
FindGroup := function(a, p)
  // Find the primes we're interested in (and -1, which we call a prime).
  S := Setseq( &join { Seqset(PrimeDivisors(a[i])) : i in {1..3} } 
                join { -1, 2 } );

  // Now form \sqrt[4]{Q_S} / Q_S.  This has one generator of order 4
  // for each element of S, including -1.
  Q4 := AbelianGroup([ 4 : p in S ]);

  // If p <> 0, compute the group \sqrt[4]{Q_p} / Q_p* and the map
  // from Q4 to it
  f := LocalFourthRoots(p, Q4, S);

  // Create Y
  Y := AbelianGroup([4,4,4,2,4]);

  // Create the evaluation map b
  b := hom<Y->Q4 | 
         [ Q4 ! GetRoot(x,S,4) : x in [ a[i] : i in {1..3} ]] 
         cat [ Q4 ! GetRoot(2,S,2), Q4 ! GetRoot(-4,S,4) ] > * f;

  // Let Z be its codomain
  // (either \sqrt[4]{Q_S}/Q_S or \sqrt[4]{Q_p}/Q_p* )
  Z := Codomain(b);

  // Compute the duals of all the groups...
  C4 := CyclicGroup(GrpAb, 4);
  YDual, tY, uY := Hom2(Y, C4);
  ZDual, tZ, uZ := Hom2(Z, C4);

  // ... and the map b
  bdual := hom< ZDual->YDual | 
    [ uY(b * tZ(ZDual.i)) : i in [1..Ngens(ZDual)] ] >; 

  // We also look at the subgroup generated by Y.5 (aka 1+i)
  // and the dual of the inclusion map.  The kernel of the dual
  // of the inclusion map consists of those homomorphisms taking
  // 1+i to 1, so it's the image of the Abelian subgroup G' < G.
  YYDual, tYY, uYY := Hom2(sub<Y | Y.5>, C4);
  idual := hom< YDual->YYDual |
    [ uYY(tY(YDual.i)) : i in [1..Ngens(YDual)] ]>;

  // Find the image of b*
  HH := Image(bdual);

  // Construct a left inverse to phi
  v := [1,1,1,2,1];
  p := map<YDual -> G | 
    h :-> &* [ G.(6-i) ^ (Eltseq(Y.i @ tY(h))[1] div v[i]) : i in [1..5] ] >;

  // Find an element of H with h(Y.5) = -i if there is one ...
  Z := { h : h in Generators(HH) | Y.5 @ tY(h) eq 3*C4.1 } join
        { -h : h in Generators(HH) | Y.5 @ tY(h) eq C4.1 };

  // ... and find its inverse image under phi
  g5 := IsEmpty(Z) select Identity(G) else p(Rep(Z));

  // Take inverse images of the other generators too, and return
  // the subgroup they generate.
  return sub<G | { p(h) : h in Generators(HH meet Kernel(idual)) }, g5 >;
end function;