/* Assorted code related to Galois groups of local fields */


/* Compute local galois groups, assuming pol irreducible over Q_p 
   This function is for degrees < 8.
   fullentry is the standard format of the database, assuming
   all of the easier entries are filled in.  In particular we need
   the following entries:
     3: e
     4: f
    14: disc root field
    19: twin polynomial for sextics
    22: number of automorphisms
*/
localgalois(pol, p, fullentry) =
{
  local(cnt,iseven, res, resf, dgs, f2, f3);

/* Unramified, so cyclic (of degree of the field) */
  if(fullentry[3]==1,
    return(dTt(fullentry[4], 1)));
  iseven= (fullentry[14] == [1,0]);  /* use disc. root field */
/* Degrees 1 and 2 are easy enough */
  if(poldegree(pol)==1, return(dTt(1,1)));
  if(poldegree(pol)==2, return(dTt(2,1)));
/* Degree 3, just check if it is even */
  if(poldegree(pol)==3,
    if(fullentry[22]>1,
      return(dTt(3,1)),
      return(dTt(3,2))));
/* Degree 4, use automorphisms and parity */
  if(poldegree(pol)==4,
    cnt=fullentry[22];
    if(cnt==2, return(dTt(4,3)));
    if(cnt==1, if(iseven==0, return(dTt(4,5)), return(dTt(4,4))),
      if(iseven==0, return(dTt(4,1)), return(dTt(4,2)))));
/* Degree 5, automorphisms and parity suffice */
  if(poldegree(pol)==5,
    if(fullentry[22]>1,
      return(dTt(5,1)));
    if(iseven==0, return(dTt(5,3)), return(dTt(5,2))));
/* Degree 6, assumes a twin poly is in entry 19 */
  if(poldegree(pol)==6,
    res=fullentry[19];
    resf=factorpadic(res,p,50);
/* Degrees of the factors of the twin poly, which we then sort */
    dgs=vector(length(resf[,1]~), h, poldegree(resf[h,1]));
    dgs=vecsort(dgs);
    if(length(dgs) == 1, /* Twin is irreducible */
      if(iseven==1, return(dTt(6,10)), return(dTt(6,13))));
    if(dgs==[1,2,3],
    /* Check if the cubic factor is A_3 or S_3 */
      f2=resf[1,1];
      if(poldegree(f2)<3, f2=resf[2,1];
        if(poldegree(f2)<3,f2=resf[3,1]));
      f2=lift(f2*Mod(1,p^50));
      if(getdisc(f2,p)==[1,0], return(dTt(6,1)), return(dTt(6,3))));
    /* Two cubic factors, check both */
    if(dgs==[3,3],
      f2=lift(resf[1,1]*Mod(1,p^50));
      f3=lift(resf[2,1]*Mod(1,p^50));
      if((getdisc(f2,p)!=[1,0]) && (getdisc(f3,p)!=[1,0]),
        return(dTt(6,9)), return(dTt(6,5))));
    if(dgs==[2,4], if(iseven, return(dTt(6,7)));
      f2=resf[1,1];
      if(poldegree(f2)==2, f2=resf[2,1]);
      f2=lift(f2*Mod(1,p^50));
      if(getdisc(f2,p)==[1,0], return(dTt(6,6)), return(dTt(6,11))));
    if(dgs==[1,1,4], if(iseven, return(dTt(6,4)), return(dTt(6,8))));
    if(dgs==[1,1,1,3], return(dTt(6,2)));
    return([0,0,0]));
/*  This works for septics, but we don't use it since there are easier
    ways for septics over local fields */
  if(poldegree(pol)==7,
    if(fullentry[22]>1,
      return(dTt(7,1)));
    if(iseven==1, return(dTt(7,3)));
    res=myseptres(pol);
    resf=factorpadic(res,p,50);
    if(length(resf[,1]~)==4, return(dTt(7,2)), return(dTt(7,4))));
  return([0,0,0])
  }

/* New Galois group scheme [|G|, type, [data], print]
   print says how to print this group, some types may have a default
   The data depends on the type.
   The known types are:
       "t" for T notation of the current degree
       "i" for an intransitive group with background data
       "a" for an abstract transitive group, no t-number available
       "o" for just the order of the group
       "tu" for C_t.C_u partial info
       "e" for the trivial group
*/

/* Generate new codes for transitive subgroups of S_n.
   The point here is to case out ones where we have a preferred name,
   otherwise a generic fallback is used.
   gord is rarely supplied, but it is the order of the group
*/
dTt(n, t, gord=0) =
{
  if(n==1, return([n, "e", [], "&lt;<i>e</i>&gt;"]));
  if(t==1, return([n, "t", [n, 1], Str("<i>C</i><sub>" n "</sub>")]));
  if(n==3 && t==2, return([6, "t", [n,t], "<i>S</i><sub>3</sub>"]));
  if(n==4 && t==2, return([4, "t", [n,t], "<i>V</i><sub>4</sub>"]));
  if(n==4 && t==3, return([8, "t", [n,t], "<i>D</i><sub>4</sub>"]));
  if(n==4 && t==4, return([12, "t", [n,t], "<i>A</i><sub>4</sub>"]));
  if(n==4 && t==5, return([24, "t", [n,t], "<i>S</i><sub>4</sub>"]));
  if(n==5 && t==2, return([10, "t", [n,t], "<i>D</i><sub>5</sub>"]));
  if(n==5 && t==3, return([20, "t", [n,t], "<i>F</i><sub>5</sub>"]));
  if(n==6 && t==2, return([6, "t", [n,t], "<i>S</i><sub>3</sub>"]));
  if(n==6 && t==3, return([12, "t", [n,t], "<i>D</i><sub>6</sub>"]));
  if(n==6 && t==4, return([12, "t", [n,t], "<i>A</i><sub>4</sub>"]));
  if(n==6 && t==5, return([18, "t", [n,t], "<i>S</i><sub>3</sub><i>C</i><sub>3</sub>"]));
  if(n==6 && t==6, return([24, "t", [n,t], "<i>A</i><sub>4</sub><i>C</i><sub>2</sub>"]));
  if(n==6 && t==7, return([24, "t", [n,t], "<i>S</i><sub>4</sub><sup>+</sup>"]));
  if(n==6 && t==8, return([24, "t", [n,t], "<i>S</i><sub>4</sub>"]));
  if(n==6 && t==9, return([36, "t", [n,t], "<i>S</i><sub>3</sub><sup>2</sup>"]));
  if(n==6 && t==10, return([36, "t", [n,t], "<i>C</i><sub>3</sub><sup>2</sup> : <i>C</i><sub>4</sub>"]));
  if(n==6 && t==11, return([48, "t", [n,t], "<i>S</i><sub>4</sub><i>C</i><sub>2</sub>"]));
  if(n==6 && t==13, return([72, "t", [n,t], "<i>C</i><sub>3</sub><sup>2</sup> : <i>D</i><sub>4</sub>"]));
  if(n==7 && t==2, return([14, "t", [n,t], "<i>D</i><sub>7</sub>"]));
  if(n==7 && t==3, return([14, "t", [n,t], "<i>C</i><sub>7</sub>:<i>C</i><sub>3</sub>"]));
  if(n==7 && t==4, return([14, "t", [n,t], "<i>F</i><sub>7</sub>"]));
  if(n==8 && t==2, return([8, "t", [n,t], "<i>C</i><sub>4</sub><i>C</i><sub>2</sub>"]));
  if(n==8 && t==3, return([8, "t", [n,t], "<i>C</i><sub>2</sub><sup>3</sup>"]));
  if(n==8 && t==4, return([8, "t", [n,t], "<i>D</i><sub>4</sub>"]));
  if(n==8 && t==5, return([8, "t", [n,t], "<i>Q</i><sub>8</sub>"]));
  if(n==8 && t==6, return([16, "t", [n,t], "<i>D</i><sub>8</sub>"]));
  if(n==8 && t==9, return([16, "t", [n,t], "<i>D</i><sub>4</sub><i>C</i><sub>2</sub>"]));
  if(n==8 && t==12, return([24, "t", [n,t], "<i>SL</i><sub>2</sub>(3)"]));
  if(n==8 && t==13, return([24, "t", [n,t], "<i>A</i><sub>4</sub><i>C</i><sub>2</sub>"]));
  if(n==8 && t==14, return([24, "t", [n,t], "<i>S</i><sub>4</sub>"]));
  if(n==8 && t==23, return([48, "t", [n,t], "<i>GL</i><sub>2</sub>(3)"]));
  if(n==8 && t==24, return([48, "t", [n,t], "<i>S</i><sub>4</sub><i>C</i><sub>2</sub>"]));
  if(n==8 && t==25, return([56, "t", [n,t], "<i>C</i><sub>2</sub><sup>3</sup>:<i>C</i><sub>7</sub>"]));
  if(n==8 && t==36, return([168, "t", [n,t], "<i>C</i><sub>2</sub><sup>3</sup>:(<i>C</i><sub>7</sub>:<i>C</i><sub>3</sub>)"]));
  if(n==9 && t==2, return([8, "t", [n,t], "<i>C</i><sub>3</sub><sup>2</sup>"]));
  if(n==9 && t==3, return([18, "t", [n,t], "<i>D</i><sub>9</sub>"]));
  if(n==9 && t==4, return([18, "t", [n,t], "<i>S</i><sub>3</sub><i>C</i><sub>3</sub>"]));
  if(n==9 && t==5, return([18, "t", [n,t], "<i>C</i><sub>3</sub><sup>2</sup> : <i>C</i><sub>2</sub>"]));
  if(n==9 && t==8, return([36, "t", [n,t], "<i>S</i><sub>3</sub><sup>2</sup>"]));
  if(n==9 && t==9, return([36, "t", [n,t], "<i>C</i><sub>3</sub><sup>2</sup> : <i>C</i><sub>4</sub>"]));
  if(n==9 && t==16, return([72, "t", [n,t], "<i>C</i><sub>3</sub><sup>2</sup> : <i>D</i><sub>4</sub>"]));
  if(n==10 && t==2, return([10, "t", [n,t], "<i>D</i><sub>5</sub>"]));
  if(n==10 && t==3, return([20, "t", [n,t], "<i>D</i><sub>10</sub>"]));
  if(n==10 && t==4, return([20, "t", [n,t], "<i>F</i><sub>5</sub>"]));
  if(n==10 && t==5, return([40, "t", [n,t], "<i>F</i><sub>5</sub><i>C</i><sub>2</sub>"]));
  if(n==10 && t==6, return([50, "t", [n,t], "<i>C</i><sub>5</sub>wr<i>C</i><sub>2</sub>"]));
  if(n==11 && t==2, return([22, "t", [n,t], "<i>D</i><sub>11</sub>"]));
  if(n==11 && t==3, return([55, "t", [n,t], "<i>C</i><sub>11</sub>:<i>C</i><sub>5</sub>"]));
  if(n==11 && t==4, return([110, "t", [n,t], "<i>F</i><sub>11</sub>"]));
  if(n==12 && t==2, return([12, "t", [n,t], "<i>C</i><sub>6</sub><i>C</i><sub>2</sub>"]));
  if(n==12 && t==3, return([12, "t", [n,t], "<i>D</i><sub>6</sub>"]));
  if(n==12 && t==4, return([12, "t", [n,t], "<i>A</i><sub>4</sub>"]));
  if(n==12 && t==5, return([12, "t", [n,t], "<i>C</i><sub>3</sub> : <i>C</i><sub>4</sub>"]));
  if(n==12 && t==11, return([24, "t", [n,t], "<i>S</i><sub>3</sub><i>C</i><sub>4</sub>"]));
  if(n==12 && t==14, return([24, "t", [n,t], "<i>D</i><sub>4</sub><i>C</i><sub>3</sub>"]));
  if(n==13 && t==2, return([26, "t", [n,t], "<i>D</i><sub>13</sub>"]));
  if(n==13 && t==3, return([39, "t", [n,t], "<i>C</i><sub>13</sub>:<i>C</i><sub>3</sub>"]));
  if(n==13 && t==4, return([52, "t", [n,t], "<i>C</i><sub>13</sub>:<i>C</i><sub>4</sub>"]));
  if(n==13 && t==5, return([78, "t", [n,t], "<i>C</i><sub>13</sub>:<i>C</i><sub>6</sub>"]));
  if(n==13 && t==6, return([156, "t", [n,t], "<i>F</i><sub>13</sub>"]));

  return([gord, "t", [n,t]]);  /* Order has to be filled in elsewhere */
}

/* Companion for cases known to be coded by T-number */

gett(g) = g[3][2]

/* For dynamic field generation, quick Galois group computation.
   Either this is C_p : C_d for some d, or it is a tame extension.
   Inputs:
     lis is a pair [a, b], and the group is definitely C_a.C_b,
       (but maybe we can do better).  In tame cases, [a,b] is [t,u].
     ent is a full entry for the field.  In particular, we may
       need its subfields.
     If the group is not transitive (some inertia groups), intrans
       is set to 1.
 */

qgal(lis, ent, intrans=0) =
{
  local(smord,g,myn, subgal);
  myn = ent[5]; /* gets the degree of the original extension */
  if(lis[2]==1, /* Group is cyclic */
    if(lis[1]==1, return(dTt(1,1))); /* Trivial group */
    g=dTt(lis[1],1);
    if(intrans, g[2]="i");
    return(g));
  /* Now not obviously cyclic */
  if(myn==3,
    if(lis[1]*lis[2]==3, return(dTt(3,1)));
    return(dTt(3,2)));
  if(myn==5,
    if(lis[1]*lis[2] == 5, return(dTt(5,1)));
    if(lis[1]*lis[2] == 10, return(dTt(5,2)));
    if(lis[1]*lis[2] == 20, return(dTt(5,3))));
  if(myn==6 && ent[20] != 0,
    if(lis[1]*lis[2] == 6,
      if(length(ent[20]) == 2, return(dTt(6,1)),
      return(dTt(6,2))));
    if(lis[1]*lis[2] == 12, return(dTt(6,3)));
    if(lis[1]*lis[2] ==18, return(dTt(6,5))));
  if(myn==7,
    if(lis[1]*lis[2] == 7, return(dTt(7,1)));
    if(lis[1]*lis[2] == 14, return(dTt(7,2)));
    if(lis[1]*lis[2] == 21, return(dTt(7,3)));
    if(lis[1]*lis[2] == 42, return(dTt(7,4))));

  if(myn==8,
    if(ent[18] != 0,
      if(lis[1]*lis[2]==8,
        if(length(ent[18][1][2])==1, return(dTt(8,1)));
        if(length(ent[18][2][2])==1, return(dTt(8,5)));
        smord=0;
        for(j=1,length(ent[18][2][2]),
          smord=max(smord, subdegs[4][indexofpoly(ent[18][2][2][j], subdegs[4])][8]));
        if(smord==4, return(dTt(8,2)));
        if(smord==8, return(dTt(8,4)));
      , /* group of order 16 */
        smord=0;
        for(j=1,length(ent[18][2][2]),
          smord=max(smord, subdegs[4][indexofpoly(ent[18][2][2][j], subdegs[4])][8]));
        if(smord==4, return(dTt(8,7)));
        /* Uggh, we need to use a8res -> irred means T8, else T6 */
        smord = a8degs(ent[6], ent[1]);
        if(smord==[4,4], return(dTt(8,6)));
        if(smord==[8], return(dTt(8,8)));
      )));

  if(myn==9,
    if(lis[1]*lis[2] == 9,
      if(ent[18] != 0,
        if(length(ent[18][1][2])==1, return(dTt(9,1)));
        if(length(ent[18][1][2])==4, return(dTt(9,2)))));

/*    if(lis[1]*lis[2] == 11, return(Str("<i>C</i><sub>11</sub>")));  */
    if(lis[1]*lis[2] == 18,
      if(lis[1]==9,
        return(dTt(9,3)),
        return(dTt(9,4))));
    if(lis[1]*lis[2] == 27, return(dTt(9,6)));
    if(lis[1]*lis[2] == 54, return(dTt(9,10))));

  if(myn==10,
    if(lis[1]*lis[2] == 50, return(dTt(10,6)));
    if(lis[1]*lis[2] == 40, return(dTt(10,5)));
    if(lis[1]*lis[2] == 20,
      if(lis[1]==10, return(dTt(10,3)));
      if(lis[1]==5, return(dTt(10,4)));
    );
    if(lis[1]*lis[2] == 10, /* look at galois group of quintic sub */
      if(ent[18] != 0,
        smord=0;
        for(j=1,length(ent[18][2][2]),
          smord=max(smord, subdegs[2][indexofpoly(ent[18][2][2][j], subdegs[2])][8]));
        if(smord==10, return(dTt(10,2)));
        if(smord==5, return(dTt(10,1)));
        );
    );

  ); /* end of degree 10 */

  if(myn==11,
    if(lis[1]*lis[2] == 11, return(dTt(11,1)));
    if(lis[1]*lis[2] == 22, return(dTt(11,2)));
    if(lis[1]*lis[2] == 55, return(dTt(11,3)));
    if(lis[1]*lis[2] == 110, return(dTt(11,4))));  /* end of degree 11 */

  if(myn==12,
    if(length(ent[18][1][2])==3, /* 3 quadratic subs */
      if(length(ent[18][4][2])==1, return(dTt(12,18)));
      smord = subdegs[3][indexofpoly(ent[18][2][2][1], subdegs[3])][7][3][2]; /* degree 3 galois T number */
      if(smord==1, return(dTt(12,2)));
      if(smord==2, return(dTt(12,3)));
      /* should not get here */
      return(1/0);
    );
    /* 1 quadratic sub */
    smord = subdegs[6][indexofpoly(ent[18][4][2][1], subdegs[6])][7][3][2]; /* degree 6 galois T number */
    if(smord == 2, return(dTt(12,5)));
    if(smord == 5, return(dTt(12,19)));
    subgal = subdegs[4][indexofpoly(ent[18][3][2][1], subdegs[4])][7][3][2]; /* degree 4 galois T number */
    if(smord == 1,
      if(subgal == 1, return(dTt(12,1)));
      if(subgal == 3, return(dTt(12,14))));
    if(smord == 3,
      if(subgal == 1, return(dTt(12,11)));
      if(subgal == 3, return(dTt(12,12)))));   /* end of degree 12 */

  if(myn==13,
    if(lis[1]*lis[2] == 13, return(dTt(13,1)));
    if(lis[1]*lis[2] == 13*2, return(dTt(13,2)));
    if(lis[1]*lis[2] == 13*3, return(dTt(13,3)));
    if(lis[1]*lis[2] == 13*4, return(dTt(13,4)));
    if(lis[1]*lis[2] == 13*6, return(dTt(13,5)));
    if(lis[1]*lis[2] == 13*12, return(dTt(13,6))));  /* end of degree 13 */

  /* Default fallback */
  return([lis[1]*lis[2], "tu", [lis[1], lis[2]], Str("<i>C</i><sub>"lis[1]"</sub>.<i>C</i><sub>"lis[2]"</sub>")]);
}




/* A few resolvents, first sextic twins */

mytwin(pol) =
{
  local(tmp, ttmp);
  tmp=pol;
  while(issquarefree(ttmp=round(twin(tmp)))==0, tmp=poltschirnhaus(pol));
  ttmp
}

from6rts(r) = r[1]^2*r[5]^2*(r[2]*r[4] + r[3]*r[6]) + \
  r[2]^2*r[4]^2*(r[1]*r[5] + r[3]*r[6]) + \
  r[3]^2*r[6]^2*(r[1]*r[5] + r[2]*r[4]) + \
  r[1]^2*r[6]^2*(r[2]*r[5] + r[3]*r[4]) + \
  r[2]^2*r[5]^2*(r[1]*r[6] + r[3]*r[4]) + \
  r[3]^2*r[4]^2*(r[1]*r[6] + r[2]*r[5]) + \
  r[1]^2*r[3]^2*(r[2]*r[6] + r[4]*r[5]) + \
  r[2]^2*r[6]^2*(r[1]*r[3] + r[4]*r[5]) + \
  r[4]^2*r[5]^2*(r[1]*r[3] + r[2]*r[6]) + \
  r[1]^2*r[4]^2*(r[2]*r[3] + r[5]*r[6]) + \
  r[2]^2*r[3]^2*(r[1]*r[4] + r[5]*r[6]) + \
  r[5]^2*r[6]^2*(r[1]*r[4] + r[2]*r[3]) + \
  r[1]^2*r[2]^2*(r[3]*r[5] + r[4]*r[6]) + \
  r[3]^2*r[5]^2*(r[1]*r[2] + r[4]*r[6]) + \
  r[4]^2*r[6]^2*(r[1]*r[2] + r[3]*r[5])


twin(pol) =
{
  local(rts, rts2, rtsc);
  rts = polroots(pol);
  rtsc = rts;
  rts2 = vector(6,h,0);
  rts2[1] = from6rts(rts);
  for(j=2,6, rts=rtsc;
    tmp=rts[1];
    rts[1]=rts[j];
    rts[j]=tmp;
    rts2[j] = from6rts(rts));
  prod(j=1,6,x-rts2[j])
}

/* We don't use it, but here is the degree 35 resolvent for septics */

septres(pol) =
{
  local(rts);
  rts = polroots(pol);
  prod(j=1,5,prod(k=j+1,6,prod(l=k+1,7, x-(rts[j]+rts[k]+rts[l]))))
}

myseptres(pol) =
{
  local(tmp, ttmp);
  tmp=pol;
  while(issquarefree(ttmp=round(septres(tmp)))==0,tmp=poltschirnhaus(pol));
  real(ttmp)
}

a8(evenpol) =
{
  local(a,b,c,d);
  a = polcoeff(evenpol,6);
  b = polcoeff(evenpol,4);
  c = polcoeff(evenpol,2);
  d = polcoeff(evenpol,0);
  a^8 - 16*a^6*b + 96*a^4*b^2 - 256*a^2*b^3 + 256*b^4 - 128*a^4*d + 
 1024*a^2*b*d - 2048*b^2*d + 4096*d^2 + 8*a^7*x - 96*a^5*b*x + 
 384*a^3*b^2*x - 512*a*b^3*x + 128*a^4*c*x - 1024*a^2*b*c*x + 2048*b^2*c*x + 
 1536*a^3*d*x - 6144*a*b*d*x + 8192*c*d*x + 28*a^6*x^2 - 240*a^4*b*x^2 + 
 576*a^2*b^2*x^2 - 256*b^3*x^2 + 512*a^3*c*x^2 - 2048*a*b*c*x^2 + 
 4096*c^2*x^2 + 1280*a^2*d*x^2 - 7168*b*d*x^2 + 56*a^5*x^3 - 320*a^3*b*x^3 + 
 384*a*b^2*x^3 + 768*a^2*c*x^3 - 1024*b*c*x^3 - 2560*a*d*x^3 + 70*a^4*x^4 - 
 240*a^2*b*x^4 + 96*b^2*x^4 + 512*a*c*x^4 - 2176*d*x^4 + 56*a^3*x^5 - 
 96*a*b*x^5 + 128*c*x^5 + 28*a^2*x^6 - 16*b*x^6 + 8*a*x^7 + x^8
}

isevenpol(poly) =
{
  for(j=1,poldegree(poly)-1,
    if((j%2 !=0) && polcoeff(poly,j) != 0, return(0)));
  return(1)
}

/* Does a8 resolvant and returns the type */

a8degs(pol, p) =
{
  local(a8pol, fp, degs);

  if(isevenpol(pol) == 0, print("Sorry, internal error in a8type.");return([1,2,3,4,5]));
  a8pol = a8(pol);
  fp=factorpadic(a8pol,p,100);
  degs = vector(length(fp[,1]), h, poldegree(fp[h,1]));
  if(sum(j=1,length(degs), degs[j]) != 8, print("Sorry, non-separable resolvent."); return([1,2,3,4,5]));
  return(vecsort(degs));
}