/* If you just want polynomials, you can comment these out and use
   quickdegp function below */

read("panayi.gp");
read("localgalois.gp");

/* Routines here are for systematically generated p-adic extensions,

   The main routine is doitall(p,n) which produces all of the data we
   generate for the database in our "standard" format.  It can be used
   when p=n or when gcd(p,n)=1.  Along the way, it generates the same
   data for proper non-trivial divisors of n to help with subfields.

   If you just want polynomials, you might use one of the following
   routines.  They can also be useful if you want say tame extensions
   of Q_2 of degree 12 withe e=3.  These do not produce full format.
   What they do produce is documented below.  See the definition of
   "doitall" near the end of the file if you want to beef their output
   up to standard format.

    unrampoljr(n,p) -> the unramified extension, just the polynomial
    tames(p,f,e)    -> the tamely ramified extensions
    padicdegp(p)    -> the ramified degree p extensions of Qp
    quickdegp(p,n)  -> puts together a couple of these routines to throw
                       in the unramified extension with the others.

   Format of the output of most of these routines is a list of
   the form [poly, c, e, f, t, slopes, d-galois, d-inertia, unram, ram]
   Here, t is tame degree of the Galois closure, slopes is a list of wild
   slopes.  The Galois group and inertia groups just give degrees (p-part
   and tame part).  Finally, it is helpful in other computations, such
   as Panayi's algorithm, to represent the extension as a totally ramified
   extension of an unramified extension.  The unramified part part is in
   the variable t and follows our standard choice of polynomial.  The
   ramified part is an Eisenstein polynomial in x or y with coefs which
   are polynomials in t.

   To convert this information to the same format used for wild extensions
   of composite degree, it takes a couple of steps.  First, use conv1 
   to do some basic conversions (syntax explained by the definition of conv1.
   To support this, code is included for discriminant root fields, and
   local root numbers.

   At the end are rountines for

     conway(n,p)   -> compute the Conway polynomial

   Routines provided here are provided without warrenty or much
   documentation.  If you want an easy interface to most of them,
   use the web front-end to the database.
*/


/* A favorite utility, we wish gp had an built in version of this
   The first input is a string giving the function to apply to each
   element of the list "lis".  The string gives the function in terms
   of z, such as mapz("z^2", [1,2,3,4,5])

   Note, you cannot nest instances of this function because the "special"
   function forbarbaz1 will have a name clash from each instance.
*/

mapz(funcstring,lis)=(forbarbaz1(z)=eval(Str(funcstring)));vector(length(lis),h,forbarbaz1(lis[h]))


/* Working toward polynomials for unramified extensions */

/* Polynomial mod p is irreducible */
/* Lots of other ways to do this, like polisirreducible(pol*Mod(1,p)) */

isirredmodp(pol, p) = poldegree(pol) == poldegree(factormod(pol,p)[1,1])

/* For an irreducible polynomial mod p: is it primitive? */
isprimpol(pol,p) =
{
  local(m,oo, pps);

  m=poldegree(pol);
  pol=pol*Mod(1,p);
  oo=p^m-1;            /* Order of the group */
  pps=factor(oo)[,1]~; /* primes dividing the order of the group */
  for(j=1,length(pps),
    if(lift(lift(Mod(x,pol)^(oo/pps[j]))) == 1, return(0)));
  return(1);
}

/* Unramified defining polynomial matching the criteria in our paper
   The only slightly tricky part is stepping through the proper lex
   ordering on polynomials. */
unrampoljr(m, p) =
{
  local(trying, pol, j,s);
  trying=1;
  pol=x^m;
  while(trying,
    j=0; s=1;
    while(polcoeff(pol,j)==(p-1)*s, pol -= s*(p-1)*x^j;s= -s;j++);
    pol += s*x^j;
    if(isirredmodp(pol, p), if(isprimpol(pol,p), trying=0)));
  pol
  }

/* Generating tame polynomials */

/* Utility for next function */

removezeros(lis) =
{
  local(cnt=0, ans);

  for(j=1,length(lis),if(lis[j],cnt++));
  ans=vector(cnt);
  cnt=0;
  for(j=1,length(lis),if(lis[j],cnt++;ans[cnt]=lis[j]));
  ans
  }
/* For tame routine, try to lift an orbit mod g to one mod p^f-1
   of length f */

liftorb(orb,p,f,g) =
{
  local(m, v,s,t);
  m=p^f-1;
  v=vector(f,h,0);
  for(k=0,m/g-1,
    for(j=1,length(orb),
      s=orb[j]+k*g;
      v[1]=s;
      t=1;s=(v[1]*p) % m;
      while(v[1] != s, t++;v[t]=s;s=(v[t]*p) % m);
      if(t==f, return(v))));
  return(-1);
  }

/* Produce a list of tame extensions with the requested data */
tames(p,f,e) =
{
  local(unrp, ans, g,v,vorb,k,norb, f2, g2, fm);

  if(e==1, unrp=unrampoljr(f,p);
    return([[unrp, 0,1,f,1,[],[f,1],[1,1],subst(unrp,x,t),y-p]]));
  unrp = subst(unrampoljr(f,p),x,y);
  g = gcd(e, p^f-1);
  /* Extra degree needed to get roots of 1 */
  f2 = lcm(f,znorder(Mod(p,e)))/f;   
  /* Compute needed orbits */
  v=vector(g-1,h,1);
  vorb=[[0]];
  for(j=1,g-1,
    if(v[j]!=0,
      norb=[j];
      k=(j*p) % g; v[j]=0;
      while(k!=j, v[k]=0; norb=concat(norb,[k]); k=(k*p) % g);
      vorb=concat(vorb,[norb]);
      ));
  /* print("f = ", f, " g = ", g, " Orbits ", vorb); */
  /* Now look for lifts */
  v=vector(length(vorb),h,liftorb(vorb[h],p,f,g));
  /* print("Lifted to: ", v); */
  /* Strategy one */
  for(j=1,length(v),
    /* Compute extra unramified part for this one */
    fm = 1;
    gm = gcd(p^(f*f2*fm)-1,e);
    while((vorb[j][1]*(p-1)*(p^(f*f2*fm)-1)/(p^f-1) % gm) != 0,
      fm++;
      gm = gcd(p^(f*f2*fm)-1,e));
    
    /* Now compute the polynomial */
    if(v[j]!= -1,
      norb = polresultant(x^e-p*y^v[j][1], unrp,y);
      norb /= pollead(norb);
      v[j] = [norb, f*f2*fm, v[j][1]],
      /* Root shift procedure */
      norb = polresultant((x+y)^e-p*y^vorb[j][1],unrp,y);
      norb /= pollead(norb); /* be safe */
      v[j] = [norb, f*f2*fm,vorb[j][1]];
      ));
  for(j=1,length(v),
    v[j]=  [v[j][1],f*(e-1),e,f,e,[],[e, v[j][2]],[e,1],
                subst(unrp,y,t),y^e-p*t^v[j][3]]);
  return(v);
  }

/* Generate all degree p ramified extensions of Qp */
/* Format is poly, c, e, f, t, slopes, d-galois, d-inert */
padicdegp(p) =
{
  local(ans, c, d1, d2, g,m, cnt,unr);

  cnt=0;
  ans=vector(p^2,h,0);
  unr = subst(unrampoljr(1,p),x,t);
  for(j=0,p-1,
    ans[cnt++]=[x^p+p*(1+j*p), 2*p-1, p, 1, p-1, [(2*p-1)/(p-1)], 
                    [p, p-1], [p,p-1], unr,y^p+p*(1+j*p)]);
  for(j=0,p-1,
    ans[cnt++]=[x^p-p*x^(p-1)+p*(1+j*p),2*p-2,p,1,1,[2],[p,1],
                     [p,1],unr,y^p-p*y^(p-1)+p*(1+j*p)]);
  for(j=1,p-2,
    c = 2*p-2;
    g=gcd(p-1,c);
    m = znorder(Mod(-j, p));
    d2= (p-1)/gcd((p-1)/m, g);
    ans[cnt++]=[x^p+j*p*x^(p-1)+p,2*p-2,p,1,1,[2],[p,d2],
                      [p,1],unr,y^p+j*p*y^(p-1)+p]);
  for(j=1,p-1,
    for(k=1,p-2,
      c = p+k-1;
      g=gcd(p-1,c);
      m = znorder(Mod(k*j, p));
      d1=(p-1)/g;
      d2= (p-1)/gcd((p-1)/m, g);
    ans[cnt++]=[x^p+j*p*x^k+p,c,p,1, d1, [c/(p-1)],
                     [p,d2],[p,d1],unr,y^p+j*p*y^k+p]));
  ans
  }

/* Put these together for easier access */
quickdegree(p,n) =
{
 local(pls=[],unrp);
  if(p==n,
    unrp = unrampoljr(p,p);
    pls = concat([[unrp,0,1,p,1,[],[p,1],[1,1],subst(unrp,x,t),y-p]],
      padicdegp(p)),
    dn=divisors(n);
    for(j=1,length(dn), pls = concat(pls, tames(p,dn[j],n/dn[j])));
  );
  pls
}

/* Start of conversion routines from the simple format used for dynamic
   entries to the format used for wild extensions of composite degree.
   Inputs are the single entry, and the prime p.  So if you have
      p=5; n=6
      mylist = quickdegree(p,n)
   Then you would need
      mylist2 = mapz("conv1(z, p)", mylist)
   to get the list.
*/

/* For easy data */

conv1(ent, p) =
{
  local(ent2, pdc,uu);
  pdc=pdiscclass(poldisc(ent[1]), p);
  uu = ent[7][1]*ent[7][2]/ent[8][1]/ent[8][2];
  ent2= [p, ent[2], ent[3], ent[4], ent[3]*ent[4], ent[1], 0, 0,
    0, ent[6], [ent[5], uu] , [ent[9],ent[10]], 0, pdc, pdc[1], pdc[2],
    dohwp(ent[1],p,ent[3]), 0, 0, 0, -1, 0,ent[7], ent[8]];
  /* ent[7] and ent[8] are in 23 and 24 */
  /* special normalization for sorting on 1st part of disc field */
  if(p==2 && ent2[15]<0, ent2[15] *= -3/2);
  ent2[21] = slopeave(p, ent[5], ent[6]); /* compute GMS */
  /* ent2[22] = countroots(ent2[6], ent2[6], p, ent2[12]);  automorphisms */
  return(ent2);
}

/* Take a discriminant and a prime p and give the local discriminant
   root field in our internal format */

pdiscclass(disc,p) =
{
  local(tdis, expo);
  tdis = disc;
  expo=valuation(tdis,p); tdis = tdis/p^expo;
  expo = expo % 2;
  if(p==2, tdis = tdis % 8, tdis = tdis % p);
  if(p==2, if(tdis==1, [p^expo, 0], if(tdis==5, [p^expo, 1], 
     if(tdis == 7, [- p^expo, 0], [- p^expo, 1]))), 
  [p^expo, if(kronecker(tdis,p)==1,0,1)])
}

/* Wrapper for easier access to pdiscclass from a polynomial and the
   prime p */

getdisc(poly, p) = pdiscclass(poldisc(poly), p)

/* Local root numbers  */

/* Hasse-Witt invariant for just for one prime p.
   Deal with easy cases, otherwise use orthogonalization.

   This only gives the adjustment from the invariant for the
   quadratic field.  See below for the final part.
*/

dohwp(inpol, p, e) =
{
  local(n, fudge, delt, f, d);
  if((p % 2)==0, 
    /* Odd degree is easy */
    if((poldegree(inpol) % 2) ==1, return((-1)^((e^2-1)/8)),
    return(dohw(inpol, [p])[1])));
  /* Now p must be odd */
  n=poldegree(inpol);
  delt = (-1)^(n*(n-1)/2);
  d = poldisc(inpol);
  fudge = hilbert(delt, d, p);
  if((n % 2)==1, return(fudge));
  if((e % p) != 0, /* Tame */
    if((e % 2)==1, return(fudge));
    f = n/e;
    if((p % 4) == 1,
      if((f % 2) == 0, return(-hilbert(p, d, p)));
        return(hilbert(p, e, p))); /* f odd */
    /* Now p === 3 (mod 4) */
      if((f % 2) == 1, return((-1)^(f*(f+1)/2)*hilbert(p, e, p)*hilbert(d,delt,p)));
        return((-1)^(f*(f+1)/2)*hilbert(p*d, p*delt, p)); /* f even */
  );
  /* Fall back */
  return(dohw(inpol, [p])[1])
}

/* dohw(poly, plist) gives a vector of invariants for the list of
   primes (code originally written for number fields.
   0 represents the infinite place.  If you want only one
   prime p, then dohw(poly, [p])[1] does the trick.
*/

dohw(inpol, rplist)=inv2(ma(inpol), rplist)

raml(list, rplist) =
{
  vector(length(rplist), h, prod(k=1, length(list)-1, 
    prod(j=k+1, length(list), hilbert(list[k],list[j], rplist[h]))))
}

inv2(m, rplist) = raml(gs(m), rplist)

ma(inpol) = matrix(length(inpol)-1,length(inpol)-1,xx,yy, trace(Mod(x,inpol)^(xx+yy-2)))

gs(m) = 
{
  local(len, ttmp, gotit);
  gotit = gramschmidt(m);
  if(gotit == 0, len=length(m);
    ttmp = matrix(len,len,xx,yy, (random() % 11)-5);
  if(matdet(ttmp), gs(ttmp*m*(ttmp~)), gs(m)), gotit)
}

gramschmidt(m,u,w,ttmp, ok) = 
{
  u = vector(length(m),j,vector(length(m),k, j==k));
  w=u; ok=1;
  for(j=1, length(m), 
    if(j>1 && w[j-1]*m*w[j-1]~ == 0, return(0));
      w[j] = u[j]-sum(k=1,j-1,(u[j]*m*w[k]~)/(w[k]*m*w[k]~) *w[k]);
  );
  w = matrix(length(m),length(m),xx,yy,w[xx][yy]);
  ttmp = w*m*(w~);
  vector(length(m), j, ttmp[j,j]*denominator(ttmp[j,j])^2)
}

/* To take the local root number from the form it is stored and print
   it, we use prinhw.  Supporting routines follow it.
   Here, p is the prime, invec is the discrim. root field (in the
   internal format as a pair), and j is the adjustment value computed
   by dohw above.  If we have an entry "ent", then the call is
     prinhw(ent[1], ent[14], ent[17])
 */

prinhw(p,invec,j) = 
{
  local(temp);
  temp = (1-hilbert(2,myquaddisc(invec,p),p)+ qhw(invec,p)+(1-j)) % 4;
  prinhwsimple(temp);
}

/* Converts a quadratic discriminant root field our from internal format
   to an integer d s.t. the field is Q_p(sqrt(d)) */
myquaddisc(invec,p) =
{
  local(temp);
  if(invec[2]==0, 
    quaddisc(invec[1]), 
     /* unramified twist */
    if(p==2, return(quaddisc(5*invec[1])));
    temp=1; 
    while(kronecker(prime(temp),p) != -1, temp++);
    quaddisc(prime(temp)*invec[1]))
}

/* Local root numbers for quadratic fields 
   The quadratic field comes in the internal format.  We also need to know p.
   The result is an exponent mod 4 so that it is really (-i)^(this number).
*/
qhw(tdis, p) =
{
  local(myval);

  if(tdis[1]==1, 0,
    if(p<7,
      if(p==2,
        if(tdis==[2,0], 0,
          if((tdis==[-2,0]) || (tdis==[-1,0]) || (tdis==[-1,1]), 3,
            if(tdis==[2,1], 2, 1))),
        if(p==3, if(tdis==[3,0], 1, 3),
          if(tdis==[5,1], 2, 0))),
      if((p % 4) == 3,
        if(tdis==[p,1], 3, (4-qhw(getdisc(x^2-p,2),2)) % 4),
        if(tdis[2]==0, 0, myval=1;
          while(kronecker(prime(myval), p) != -1, myval = myval+1);
          (4-(1-kronecker(-1,prime(myval)))/2- 
            qhw(getdisc(x^2-kronecker(-1,prime(myval)) 
          *prime(myval)*p,prime(myval)),prime(myval))) %4))))
  }


/* Prints (-i)^temp where i is gp's I, assuming temp is already
   reduced mod 4. */
prinhwsimple(temp) =
{
  if(temp==0, print1("1"), 
    if(temp==1, print1("-<i>i</i>"), 
      if(temp==2, print1("-1"), 
        print1("<i>i</i>"))));
}

/* Galois mean slope from the prime, the tame degree, and the list
   of wild slopes */

slopeave(p, tame, wlist) = ((p-1) * sum(j=1,length(wlist),wlist[j]*p^(j-1))+\
    (tame-1)/tame)/p^length(wlist)

/* Utility to search through subfields for field whose nearly full
   entry is given by ent, and we look through list of candidates
   lis */
findsubs(ent, lis) =
{
  local(ans=[]);

  for(j=1,length(lis),
    if((ent[3] %lis[j][3])==0 && (ent[4] %lis[j][4])==0,
      if(countroots(lis[j][6], ent[6], ent[1], ent[12], ent[13],1)>0,
        ans=concat(ans, [lis[j][6]]))));
  return(ans);
}

/* Determines degrees to check in, then finds the subs, and puts
   them in the right format */
findallsubs(ent) =
{
  local(divs);
  divs = divisors(ent[5]);
  vector(length(divs)-2,h,[divs[h+1], findsubs(ent, subdegs[divs[h+1]])])
}

/* Look up a polynomial in a pre-generated list.
   We don't try to compute any invariants of poly because it might be
   really horrible (e.g., coming from high precision p-adic factorization),
   so we take the safer route and just compare to everything.

   Returns the matching poly, or 0 if it can't find one.
*/

findinlist(poly, lis) =
{
  for(j=1,length(lis),
    if(countroots(poly, lis[j][6], lis[j][1], lis[j][12], lis[j][13],1)>0,
      return(lis[j][6])));
  return(0)
}

/* Utility, if we know a polynomial is in a list, and we want its index */

indexofpoly(poly, lis) = for(j=1,length(lis),if(lis[j][6]==poly,return(j))); 0


global(subdegs);
/* Returns completed list with full format for tame and (p,p) extensions.
   Optional parameter is so that we can call recursively to compute
   subfield data, but we don't have to keep recomputing it.
 */
doitall(p,n, havesubdegs= 0) =
{
  local(sd, mylist, fpf, twinmatches, ent);
  sd = divisors(n);
  if(length(sd)>2 && havesubdegs==0,
    /* don't need n, this has entries we won't fill, but this is easy
       to access for all degrees dividing n */
    subdegs= vector(n-1,h,[]);
    for(j=2, length(sd)-1,
      subdegs[sd[j]] = doitall(p, sd[j], 1);
    );
  );
  mylist = quickdegree(p,n);
  mylist = mapz("conv1(z, p)", mylist); /* Fast stuff */
  /* Fill in automorphisms */
  for(j=1,length(mylist),
      mylist[j][22] = countroots(mylist[j][6], mylist[j][6], p, mylist[j][12]));
  /* Fill in subfields */
  for(j=1,length(mylist),
      mylist[j][18] = findallsubs(mylist[j]));
  /* If the degree is 6, we need to compute, factor, and identify
     factors of twin.  Here, p > 5, so we only have to worry about factors
     of degree <4, which means they divide 6, which means they are already
     computed in subdegs. */
  if(n==6,
    default(realprecision, 500);
    for(j=1,length(mylist),
      mylist[j][19] = real(mytwin(mylist[j][6]));
      fpf = factorpadic(mylist[j][19], p, 500)[,1]~;
      fpf = vecsort(mapz("[z, poldegree(z)]", fpf), [2]);
      fpf = mapz("lift(Mod(1,p^500)*z[1])", fpf);
      twinmatches=[];
      for(jj=1,length(fpf),
        if(poldegree(fpf[jj])>1, /* skip degree 1 factors */
          twinmatches=concat(twinmatches,
            findinlist(fpf[jj], subdegs[poldegree(fpf[jj])]));
          ));
      mylist[j][20] = twinmatches;
    );
  );
  /* Now fill in Galois and inertia group entries.  For tame and (p,p),
     we can use quickgalois fuction (from localgalois.gp).
     Some of the data was specially stored in entries 23 and 24 by conv1 */
  for(j=1,length(mylist),
    ent = mylist[j];
    mylist[j][9] = qgal(ent[24], ent, ent[4]>1);
    mylist[j][7] = qgal(ent[23], ent);
    mylist[j][8] = ent[23][1]*ent[23][2]);
  mylist = vecsort(mylist, [2,3,15,16,21,8,22]);
  return(mylist);
}



/**********************************************************************/
/**********************************************************************/

/* For fun, Conway polynomials - but not used by db */


/* We cache ones which have been computed so far in conwaylist */
global(conwaylist);
/* A trick to initialize the cache, but only if we have not loaded this
   file already - helpful when debugging */

if(type(conwaylist)=="t_POL", conwaylist=[]);

subismult(pow, starter, modu) =
{
  local(ans);
  ans=Mod(0,modu);
  for(j=0, poldegree(starter),
    if(lift(polcoeff(starter, j))!= 0,
      ans += polcoeff(starter, j)*Mod(x^j, modu)^pow));
  lift(lift(ans)) == 0
  }

conway(m, p) =
{
  local(trying, pol, j, divs,subcon,try2);
  for(j=1,length(conwaylist),
    if(conwaylist[j][1] == [m,p], return(conwaylist[j][2])));
  trying=1;
  /* Degree 1 is basically primitive roots */

  if(m==1, return(Mod(1,p)*x-znprimroot(p)));
  
  pol=x^m* Mod(1,p);
  divs = factor(m)[,1]~;
  subcon = vector(length(divs), h, conway(m/divs[h],p)*Mod(1,p));
  while(trying,
    j=0;
    while(lift((-1)^(m-j+1)*polcoeff(pol,j))==1,
      pol -= (-1)^(m-j)*(p-1)*x^j;j++);
    if(j==m, return("failed"));
    pol += (-1)^(m-j)*x^j;
    if(isirredmodp(pol, p), if(isprimpol(pol,p),
        try2=1;
        while(0<try2 && try2<= length(divs),
          if(subismult(((p^m-1)/(p^(m/divs[try2])-1)), subcon[try2], pol),
            try2++,
            try2=-1)
          );
        if(try2>0, trying=0))));
  conwaylist = concat(conwaylist, [[[m,p], pol]]);
  pol
  }

/* End of Conway poly section */