/* cyclic_subfields.m
 * ==================
 *
 * Compute quadratic and cyclic cubic subfields of a given field.
 *
 * Author: Juergen Klueners, Stephan Elsenhans
 *
 * Last Change: May 2025
 *
 * ---------------------------------------------------------------------------- */

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

pzl := PrimesInInterval(1,1000);  
rings := [PolynomialRing(GF(p)) : p in pzl];  // This seems to be a magic speed up!

// Compute the roots of a rational polynomial in a number field by Trager's method.
function roots_by_trager(nf,pol:OnlyOne := false)
    shift := Evaluate(pol,Polynomial([-nf.1,1]));
    //norm:=Norm(shift);
      norm:=SumPoly(DefiningPolynomial(nf), pol);
//    norm, norm2;
//    assert norm eq norm2;
 fac := Factorization(norm);
 if Degree(fac[1][1]) eq Degree(pol) * Degree(nf) then
 // assert IsIrreducible(shift);
  return [];
 end if;
 if OnlyOne then
  gcd := [Gcd(fac[1][1],shift)];
 else
  gcd := [Gcd(a[1],shift) : a in fac];
 end if;
 assert &and [Degree(a) eq 1 : a in gcd];
 res := [(Roots(a)[1][1] -nf.1): a in gcd];
// assert Evaluate(pol,res[1]) eq 0; 
 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;
 assert2 IsIrreducible(f);
 assert IsPrime(deg);

 if deg eq 2 then
  r := PowerSeriesRing(Rationals(),Degree(f) div 2 + 1);
  g := Sqrt(r!Reverse(f));
//  g;
  f2 := Polynomial(Reverse([Coefficient(g,i) : i in [0..Degree(f) div 2]]));
 else 
  r := PowerSeriesRing(Rationals(),Degree(f) div deg + 1);
  g := (r!Reverse(f))^(1/deg);
  f2 := Polynomial(Reverse([Coefficient(g,i) : i in [0..Degree(f) div deg]]));  
 end if;
// f; f2; f - f2^deg;
 dd := GCD([Numerator(a) : a in Coefficients(f - f2^deg)]);
 if verb then printf"Factoring %o\n",dd; end if;
 fac := Factorization(dd);
 pzl := [a[1] : a in fac];
 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;
 last_rank_change := 0;
 while (Degree(f) mod 2^(#l - r) ne 0) or (#mat lt 3) or (last_rank_change lt 7) or (1 in rs and r lt #l) do  
  if (not p in l) then
   fp := PolynomialRing(GF(p))!f;
   
   if not IsSquarefree(fp) then p := NextPrime(p+1); continue; end if;

   last_rank_change :=  last_rank_change + 1;
   prime_count := prime_count + 1;

   fac := DistinctDegreeFactorization(fp);
   dl := [t[1] : jj in [1..Degree(t[2]) div t[1]], t in fac];   

   deg_vec := 1;
   for d in dl do
     deg_vec := BitwiseOr(deg_vec, ShiftLeft(deg_vec,d));
   end for;
   no_split := BitwiseAnd(deg_vec, ShiftLeft(1,Degree(f) div 2)) eq 0;
   if 1 in [a mod 2 : a in dl] 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(f,verb);
 if Signature(nf) 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;
  if verb then
   printf"Testing %o\n",d;
  end if;
  r := roots_by_trager(nf, x^2-d:OnlyOne);
  if verb then
   printf"%o roots found \n",r;
  end if;
  if #r gt 0 then
      //Append(~res,sub<nf | r[1][1]>);
   if #suc eq 0 then 
       suc := [d];
       kk:=NumberField(x^2-d);
       Embed(kk,nf, r[1]);
       Append(~res,kk);
   else
    n := #suc;
    for j := 1 to n do
     d1 := d * suc[j] div GCD(d,suc[j])^2;
     Append(~suc,d1);
     kk:=NumberField(x^2-d1);
     Embed(kk,nf, r[1]*res[j].1 / GCD(d,suc[j]));
     Append(~res,kk);
     //Append(~res,sub<nf | r[1][1]*res[j].1 / GCD(d,suc[j])>); 
    end for;
    Append(~suc,d);
    kk:=NumberField(x^2-d);
    Embed(kk,nf, r[1]);
    Append(~res,kk);
   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 if a system of invariant blocks of size size is compatible with the factorization pattern dl
function has_invariant_blocks(dl, size)
 Sort(~dl);
 if dl[#dl] gt size then return false; end if;
 i := 0;
 while (i lt #dl) and (dl[i+1] lt size) do i +:= 1; end while;
 if i le 2 then return true; end if;
 mm := dl[i]; 
 if (mm le 2) then
  if (size mod mm eq 0) then return true; end if;
// mm is 2 and the block size is odd
  bl_cnt := &+(dl[1..i]) div size;
  return #[a : a in dl | a eq 1] ge bl_cnt; // We need at least one 1 for each block  
 end if; 
 if dl[1] + dl[i] gt size then return false; end if;
 
 if dl[1] + dl[i] eq size then 
  dl := dl[2..i-1]; 
 else
  dl := dl[1..i];
 end if;
 
 bl_cnt := &+dl div size; // We have to build this many blocks

 if dl[1] + dl[#dl] gt size then 
  return false;
 end if;
 
 if bl_cnt le 1 then return true; end if;

 ds := {0};
 for a in dl do ds := ds join {a + t : t in ds}; end for;
 if not size in ds then return false; end if;
// We rarely get here...
 if bl_cnt eq 2 then return true; end if;
 if not 2*size in ds then return false; end if;
 
// Remaining case: We have to build 3 blocks, all cycles are smaller than the block size.
 
// We test, if each cycle can be extended to an invariant block. 
 m := Set(dl);
 for a in m do
  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;
   if size - a in sums then continue a; end if;
  end for;
  if not size - a in sums then
   return false; 
  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];
 cubs := [nf3!1];
 for p in ram do
  fac := Factorization(p*ord_3);
  assert #fac eq 2;
  e1,e2 := IsPrincipal(fac[1][1]);
  pi2 := p / e2;
  Append(~cubs,pi2^3);
  Append(~rads,e2 * pi2^2);
 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 p in ram then continue; end if;
  fp := PolynomialRing(GF(p))!f;
  if not IsSquarefree(fp) then continue; end if;
  fac := DistinctDegreeFactorization(fp);
  if &and [a[1] mod 3 eq 0 : a in fac] then continue; end if;
  dl := [a[1] : b in [1..Degree(a[2]) div a[1]],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]);  
  nk := #Kernel(Transpose(Matrix(GF(3),mat)));
 until assigned(nk) and ((#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 := [];  
 co_l := [];
// Convert the kernel to a list of radicands
 gew := [];
 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]] / &*([1] cat [cubs[i] : i in [1..#rads] | co[i] eq 2]));
   Append(~co_l,a);
   Append(~gew,&+co);
  end if;
 end for;
 Sort(~gew,~perm); // Try to keep the twists small by this type of sorting 
 res := PermuteSequence(res,perm);
 co_l := PermuteSequence(co_l,perm);
 return res, co_l;
end function;

// Convert radicands to polynomials defining cubic extensions of Q
function polys_from_radicands(res, co_l)
 pol_l := [];
 co_res := [];
 for i := 1 to #res do
  alpha := res[i];
  suc, na3 := IsPower(Norm(alpha),3);
  if not suc then continue i; end if;
// 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 i; end if; // No cubic subfield possible
   if Degree(p1) eq 3 then // Candidate of cubic field found
     Append(~pol_l,p1); 
     Append(~co_res,co_l[i]);
     continue i; 
   end if;
  end for;
 end for;
 return co_res, pol_l;
end function;


/* Construct all the twists of sf with the fields in res */
/* All the roots of the defining polynomials are stored in rt_ll, r_sf */
list3:=[1..3];
function cubic_twists(sf_l, sf_new, rt_ll, r_sf, nf)
 res := [];
 for i := 1 to #sf_l do
  gen1 := &+[(rt_ll[i][j]) * (nf!r_sf[j]) : j in list3];
  gen2 := &+[(rt_ll[i][j]) * (nf!r_sf[4-j]) : j in list3];
//  MinimalPolynomial(gen1); MinimalPolynomial(gen2); // This verbose printing is expensive
  
  for g in [gen1,gen2] do
   g2 := g*g; g3 := g2*g; // As we expect a cubic minimal polynomial, this is the fastest way to get it.
   ker := Kernel(Matrix([Eltseq(g3),Eltseq(g2),Eltseq(g),Eltseq(nf!1)]));
   assert Dimension(ker) gt 0;
   if Dimension(ker) ne 1 then
    printf"Tschirni to implemented\n";
    assert false;
   end if;
   bm := BasisMatrix(ker);
   mp := Polynomial(Reverse(Eltseq(bm)));
   assert2 Evaluate(mp,g) eq 0;
   sf := NumberField(mp);
   Embed(sf,nf,g);
   Append(~res,sf);
  end for; 
 end for;
 return res;
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 eq 1]; // Tame ramification of a cyclic extension requires p = 1 mod 3
  
 rad_1,co_l := potential_radicands(DefiningPolynomial(nf),ram);
 co_l,p3l := polys_from_radicands(rad_1, co_l);
 res := [];
 rt_ll := [];
 co_done := {};
 for i := 1 to #p3l do
  if co_l[i] in co_done then continue i; end if; // Field already found by twisting
  f3 := p3l[i];  
  if IsSquare(Discriminant(f3)) then
   r := roots_by_trager(nf,f3:OnlyOne); 
   if #r gt 0 then
    sf := NumberField(f3);
    Embed(sf,nf,r[1]); 
    r_sf := [nf!nn[1] : nn in Roots(f3,sf)];
    if #res gt 0 then
     l := cubic_twists(res, sf, rt_ll, r_sf, nf);
     for a in l do
      Append(~res,a);
      Append(~rt_ll,[nf!nn[1] : nn in Roots(DefiningPolynomial(a),a)]);
     end for;
     co_done := co_done join {c1*a + c2*co_l[i] : c1,c2 in [1,2], a in co_done}; 
    end if;
    Append(~res,sf); 
    Append(~rt_ll,r_sf);
    Include(~co_done, co_l[i]);
   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);

c3l := [
    t^3 - 7*t^2 + 14*t - 7,
    t^3 - 13*t^2 + 26*t - 13,
    t^3 - 19*t^2 + 38*t - 19,
    t^3 - 93*t^2 + 899*t - 31];
nf_ab := SplittingField(c3l);
nf_ab := NumberField(PolynomialRing(Integers())!DefiningPolynomial(nf_ab));    
ll := cyclic_cubic_subfields(nf_ab);

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

-----------

Test with database:

r<x> := PolynomialRing(Integers());
load"polys16";
#L;
res := []; 
time for f in L do 
 a := quadratic_subfields(NumberField(f)); 
 Append(~res,#a); 
end for;
[#[a : a in AllPartitions(TransitiveGroup(16,i)) | #a eq 8] : i in [1..#res]] eq res;

load"pol_18.txt";
L := pol_l;
#L;
res := []; 
time for f in L do 
 a := quadratic_subfields(NumberField(f)); 
 Append(~res,#a); 
end for;
[#[a : a in AllPartitions(TransitiveGroup(18,i)) | #a eq 9] : i in [1..#res]] eq res;

res := []; 
time for f in L do 
 a := cyclic_cubic_subfields(NumberField(f)); 
 Append(~res,#a); 
end for;
[#[a : a in AllPartitions(g) | #Image(Action(g,Orbit(g,a))) eq 3] where g is TransitiveGroup(18,i) : i in [1..#res]] eq res;

load"poly_20.txt";
#L;
res := []; 
time for f in L do 
 a := quadratic_subfields(NumberField(f));  

 Append(~res,#a); 
end for;
[#[a : a in AllPartitions(TransitiveGroup(20,i)) | #a eq 10] : i in [1..#res]] eq res;

load"poly_van_hoeij.txt";
#L;

res := []; 
time for f in L do 
 a := quadratic_subfields(NumberField(f));  

 Append(~res,#a); 
end for;
[a : a in res];

res := []; 
time for f in L do 
 a := cyclic_cubic_subfields(NumberField(f));  

 Append(~res,#a); 
end for;
[a : a in res];

*/
