/* cyclic_subfields.m
 * ==================
 *
 * Compute quadratic and cyclic cubic subfields of a given field.
 *
 * Author: Stephan Elsenhans
 *
 * Last Change: July 2019                                                    
 *
 * ---------------------------------------------------------------------------- */

/* Functions to be called:
 * =======================
 *
 * quadratic_subfields(nf) : All quadratic subfields of nf
 *
 *
 * cyclic_cubic_subfields(nf) : All cyclic cubic subfieds of nf
 *
 * ---------------------------------------------------------------------------- */



// Compute a potentiall e-th root of a polynomial
function pol_root(f,e)
 assert Degree(f) mod e eq 0;
 assert LeadingCoefficient(f) eq 1;
 n := Degree(f);
 d := n div e;
 X := PolynomialRing(Rationals()).1;
 res := X^d;
 di := f - X^Degree(f);
 for k := 1 to d do 
  term := Coefficient(di,n-k) / e * X^(d-k);
  res := res + term;
  di := f - res^e;
 end for;
 return res;
end function;

// Compute a superset of the odd ramified primes of every quadratic subfield
function potentially_ramified_primes(f,verb:deg := 2)
 assert BaseRing(Parent(f)) eq Integers();
 assert LeadingCoefficient(f) eq 1;
 assert Degree(f) mod deg eq 0;
 assert Degree(f) gt 1;
 assert IsIrreducible(f);
 assert IsPrime(deg);

 ps :=  PolynomialToPowerSums(f);
 ps2 := [ps[i] / deg : i in [1..#ps div deg]];
 f2 := PowerSumsToPolynomial(ps2);
 dd := GCD([Numerator(a) : a in Coefficients(f - f2^deg)]);
 if verb then printf"Factoring %o\n",dd; end if;
 pzl := PrimeDivisors(dd);
 if not deg in pzl then Append(~pzl,deg); end if;
 return pzl;
end function;

// Use cycle types to determine a set of potential discriminants of quadratic subfields
function cycle_test(l,f)
 l := Sort(l);
 mat := []; rs := [];
 p := NextPrime(Degree(f));
 prime_count := 0;
 r := 0;
/* while ((#mat lt #l + 5) and (prime_count lt Degree(f)*#l or #mat lt 3) and (#mat le 2* r)) or 
       (1 in rs and r lt #l) or (2^(l - r) gt Degree(f))  do */
 last_rank_change := 0;
 while (Degree(f) mod 2^(#l - r) ne 0) or (#mat lt 3) or (last_rank_change lt 10) or (1 in rs and r lt #l) do  
  if (not p in l) then
   last_rank_change :=  last_rank_change + 1;
   prime_count := prime_count + 1;
   fac := Factorization(f,GF(p));

   deg_vec := 1;
   for j := 1 to #fac do
    dd := Degree(fac[j][1]);
    for k := 1 to fac[j][2] do
     deg_vec := BitwiseOr(deg_vec, ShiftLeft(deg_vec,dd));
    end for;
   end for;
   no_split := BitwiseAnd(deg_vec, ShiftLeft(1,Degree(f) div 2)) eq 0;
   if 1 in [(Degree(a[1])*a[2]) mod 2 : a in fac] then // Prime splits in each quadratic subfield
    if no_split then
     return []; // Cycle excludes split: No quadratic subfield  
    end if;
    Append(~mat,[ (JacobiSymbol(a,p) eq -1) select 1 else 0 : a in l]);
    Append(~rs,0);
   else
    if no_split then // Cycles proves: prime is inert 
     Append(~mat,[ (JacobiSymbol(a,p) eq -1) select 1 else 0 : a in l]);
     Append(~rs,1);
    end if;
   end if;
  end if;
//  p,mat;
  p := NextPrime(p+1);
  r_old := r;
  r := Rank(Matrix(GF(2),mat)); 
  if r ne r_old then
   last_rank_change := 0;
  end if;
  if 1 in rs then
   if not IsConsistent( Transpose(Matrix(GF(2),mat)), Vector(GF(2),rs)) then return []; end if;
  end if;
 end while;
 suc, sol,ker := IsConsistent( Transpose(Matrix(GF(2),mat)), Vector(GF(2),rs));
 if not suc then return []; end if;
 if Set(rs) eq {0} then
  return [ &*([l[i] : i in [1..#l] | s[i] eq 1]) : s in ker | not IsZero(s)];
 end if;
 return [ &*([l[i] : i in [1..#l] | s[i] eq 1] cat [1]) where s is sol + a : a in ker];
end function;


// Find all quadratic subfields of a number field
function quadratic_subfields(nf:verb := false)
 if Degree(nf) mod 2 eq 1 then return []; end if;
 f := DefiningPolynomial(nf);
 pzl := potentially_ramified_primes(DefiningPolynomial(nf),verb);
 if Signature(NumberField(f)) eq 0 then
  Append(~pzl,-1);
 end if; // Otherwise all quadratic subfields are real.
 if verb then printf"Potentially ramified primes %o\n",pzl; end if;
 pos_discs := cycle_test(pzl,f);
 if verb then printf"Possible discriminants after cycle test %o\n",pos_discs; end if;
 x := Parent(f).1;
 res := [];
 reject := {};  
 suc := []; 
 for d in pos_discs do
  if (d in suc) or (d in reject) then continue; end if;
  r := Roots(x^2-d,nf);
  if #r gt 0 then
   Append(~res,sub<nf | r[1][1]>);
   if #suc eq 0 then 
    suc := [d]; 
   else
    n := #suc;
    for j := 1 to n do
     d1 := d * suc[j] div GCD(d,suc[j])^2;
     Append(~suc,d1);
     Append(~res,sub<nf | r[1][1]*res[j].1 / GCD(d,suc[j])>); 
    end for;
    Append(~suc,d);
   end if;
  else
   Include(~reject,d);
   for d1 in suc do
    Include(~reject,d * d1 div  Gcd(d,d1)^2);
   end for;
  end if;   
 end for;
 return res;
end function;

/* ------------------------  Cyclic cubic fields ----------------------- */

// Check is a system of invariant blocks is compatible with the factorization pattern dl
function has_invariant_blocks(dl, size)
 m := Set(dl);
 if Max(m) gt size then return false; end if;
 for a in m do
  if a lt size then
   i0 := Index(dl,a);
   sums := {0};
   for j := 1 to #dl do
    if j ne i0 then
     sums := sums join {dl[j] + a : a in sums};
    end if;
   end for;
   if not size - a in sums then
    return false; 
   end if;
  end if;
 end for; 

 return true;
end function;


// Cyclic cubic fields are radical extensions over Q(\zeta_3).
// This function computes a list of possible randicand. 
// The large field is given by f. Only the primes in ram can be ramified.
function potential_radicands(f,ram)
 nf3 := CyclotomicField(3);
 ord_3 := MaximalOrder(nf3);
 
 rads := [ord_3!nf3.1];
 for p in ram do
  fac := Factorization(p*ord_3);
  for a in fac do
   e1,e2 := IsPrincipal(a[1]);
   assert e1;
   Append(~rads,e2);
  end for;
 end for;

// Each prime proven to split in a cyclic cubic extension gives a linear relation for 
// the possible radicands
 mat := [];
 p := 5;
 repeat 
  p := NextPrime(p);
  if not p in ram then
   fac := Factorization(f,GF(p));
   if {(Degree(a[1])*a[2]) mod 3 : a in fac} ne {0} then
    // p is proven to split in every cyclic cubic subfield
    dl := [Degree(a[1]) : b in [1..a[2]],a in fac];
    if not has_invariant_blocks(dl,Degree(f) div 3) then // No subfields by cycle test
//     printf"No cyclic cubic field by cycle test.\n";
     return [],[];
    end if;
    fac := Factorization(p*ord_3);
    qu,pi := ResidueClassField(ord_3,fac[1][1]);
    red_rads := [pi(a)^((#qu - 1) div 3) : a in rads];
    m := [1] cat SetToSequence({a : a in red_rads | a ne 1});
    Append(~mat,[Index(m,a)-1 : a in red_rads]);  
   end if;
  end if;
  nk := #Kernel(Transpose(Matrix(GF(3),mat)));
 until (#mat ge #rads+4 and nk le Degree(f)) or (nk eq 1 and #mat gt 0);
 ker := Kernel(Transpose(Matrix(GF(3),mat)));
 if Dimension(ker) eq 0 then return [],[]; end if;
 res := [];  

// Convert the kernel to a list of radicands
 for a in ker do
  if not IsZero(a) then
   co := [Integers()!b : b in Eltseq(a)]; 
   i := 1; while co[i] eq 0 do i := i+1; end while;
   if co[i] ne 1 then continue; end if;
   Append(~res,&*[rads[i]^co[i] : i in [1..#rads]]);
  end if;
 end for;

 return res;
end function;

// Convert radicands to polynomials defining cubic extensions of Q
function polys_from_radicands(res)
 pol_l := [];
 rad := [];
 for alpha in res do
  suc, na3 := IsPower(Norm(alpha),3);
  if not suc then
// Non-cyclic cubic field possible
//   printf"Naive radical \n";
   Append(~pol_l,Polynomial([Norm(alpha),0,0,1])); 
   Append(~rad,alpha);
  else
// Treat cyclic degree 6 extension Q(\zeta_3,sqrt[3](a)) / Q 
   beta_3 := (Trace(alpha) - alpha) / alpha^2; 
   N_beta := 1/na3;
   nf2 := NumberField(Polynomial([-alpha,0,0,1]));
   gen2 := N_beta^2 * beta_3 * alpha^2 * nf2.1^2;
   assert gen2 eq  (nf2.1/na3)^2 * (Trace(alpha) - alpha);

// gen2 and nf2.1 are two conjugate elements. Derive a cubic cubic sufield via sum or product  
   for s := 0 to 4 do
    if s eq 0 then 
     p1 := MinimalPolynomial(gen2+nf2.1,Rationals());
    else
     p1 := MinimalPolynomial((gen2+s)*(nf2.1+s),Rationals());
    end if;
    if Degree(p1) gt 3 then printf"Fail in descent.\n"; continue alpha; end if; // No cubic subfield possible
    if Degree(p1) eq 3 then // Candidate of cubic field found
     Append(~pol_l,p1); 
     Append(~rad,alpha);
     continue alpha; 
    end if;
   end for;
  end if;
 end for;
 return rad, pol_l;
end function;


function cyclic_cubic_subfields(nf:verb := false)
 if Degree(nf) mod 3 ne 0 then return []; end if;
 ram := potentially_ramified_primes(DefiningPolynomial(nf),verb:deg := 3);
 ram := [p : p in ram | p mod 3 ne 2]; // A cyclic cubic extension does not ramify in primes = 2 mod 3
  
 rad_1 := potential_radicands(DefiningPolynomial(nf),ram);
 rad,p3l := polys_from_radicands(rad_1);
 res := [];
 for a in p3l do
  if IsSquare(Discriminant(a)) then
   r := Roots(a,nf);
   if #r gt 0 then
    sf := sub<nf | r[1][1]>;
    if Degree(sf) eq 3 then
     Append(~res,sf); 
    end if;
   else
    "Subfield in final state excluded";
   end if;
  end if;
 end for;
 return res;
end function;


/*  A few test examples:

r<x> := PolynomialRing(Integers());
pol16 := x^16 - 136*x^14 + 6476*x^12 - 141912*x^10 + 1513334*x^8 
          - 7453176*x^6 + 13950764*x^4 - 5596840*x^2 + 46225;

nf1 := NumberField(pol16);
time a := quadratic_subfields(nf1);    

pol81 :=  x^81 - 9*x^78 + 261*x^75 - 6465*x^72 - 5976*x^69 - 88659*x^66 + 1438305*x^63 + 142618428*x^60 
          - 1061145738*x^57 + 998018267*x^54 - 7764846651*x^51 + 16670823876*x^48 + 6321985362*x^45 
          + 308149787331*x^42 + 1399281423768*x^39 + 4617405757797*x^36 + 10048968487008*x^33 
          + 21269478032217*x^30 + 37113067712761*x^27 + 46571832382185*x^24 + 52747912349847*x^21 
          + 43747548047208*x^18 + 24276383580735*x^15 + 14013201103731*x^12 
          + 508747247625*x^9 - 147357261396*x^6 - 1879657452*x^3 - 357911;

nf2 := NumberField(pol81);
time b := cyclic_cubic_subfields(nf2);

pol54 := x^54 - 476*x^52 - 145786*x^50 + 66346872*x^48 + 12203235501*x^46 - 3042974602432*x^44 
         - 629064777910472*x^42 + 25174890008006146*x^40 + 12468641433551186616*x^38 + 1014181518140160958020*x^36 
         + 34731237414858100773108*x^34 + 453812086601387588349594*x^32 - 1490639387776131163290056*x^30 - 
    72116436587808498770054652*x^28 - 219510831493555231703735107*x^26 + 942992234071220518773275172*x^24 
    + 2849155850696579809849265810*x^22 - 1384045828336848111135015520*x^20 - 9601851091705417313624904255*x^18 
    - 9239275450683071056428277316*x^16 - 3105928150710779984291465504*x^14 - 342171584172961915759891648*x^12 
    + 213848671357767117939968*x^10 + 237631081210097994752*x^8 - 55684023253823488*x^6 - 52142535901184*x^4 
    - 7065829376*x^2 - 262144;

nf3 := NumberField(pol54);
time c := cyclic_cubic_subfields(nf3);

*/

