/* The file is devoted to Panayi's algorithm (2 implimentations).

   In almost all cases, we just need to count roots.  The main
   countroots program wants some extra optional information.
   The simplest to compute is nfinit and an ideal above p (from
   idealprimedec).  If you are comparing several polynomials to
   a fixed one, you can precompute these and feed them into the function.

   Alternatively, use entry 12 from a database field.  This gives the
   polynomial for the unramified subextension, and then an Eisenstein
   polynomial for the ramified part (in terms of the unramified part).
   In this case, the pr entry is irrelevant.  If nf is like this,
   (specifically if it has length 2), then we jump to an alternate
   version of this program.

   Finally, if you want just a yes/no decision from countroots, set
   the optional variable early to 1.
*/
/* This sets the precedence of variables for gp.  If you have done
   anything with these variables in gp before loading this file, you
   may be sunk and countroots will either spin until it runs out of
   memory, or cause an error.
*/
x;
y;
t;

/* Counting roots of f in Q_p/(g)
   We can optionally add data about g:
     nf is its number field, and pr is the prime above p we want
   This version uses nfinit, either directly, or you pass it the
   nf structure.
*/
countroots(f, g, p, nf, pr, early=0) =
{
  local(clis, rlist, cspot, rspot,m, c, gg, hh,ldeg, unif, resg, rrs);

  /* If we have new type data, use it instead */
  if(length(nf)==2, return(countrootsdb(p,f,nf[1], nf[2],early)));
  clis = vector(poldegree(f),h,0);
  rlis = vector(poldegree(f),h,-1);
  cspot=1;
  gg=subst(g,x,y);
  if(nf==0, nf=nfinit(g));
  if(pr==0, pr=idealprimedec(nf,p)[1]);
  if(nfeltval(nf, pr[1] ,pr)==1, unif=Mod(pr[1],gg),
    if(nfeltval(nf, pr[2], pr)==1, unif=Mod(subst(nf.zk*pr[2],x,y),gg),
      print("Cannot find a local uniformizer"); 1/(2-2)));
  clis[1]= psharp(f*Mod(1, gg), p, nf, pr, unif);
  m=0;
  /* Find a generator for the residue field */
  if(p==2 || pr[4]==1,
    resg=nfprimroot(nf,pr);
    );
  while(cspot>0,
    c=clis[cspot];
    cspot--;
    /* This is where we find the roots of c modulo p.
    */
    if(p==2 || pr[4]==1,
      cbar= vector(poldegree(c)+1,h,subst(lift(polcoeff(c,h-1)),y,x));
      /* Test to see if 0 is a root mod pr */
      rspot=0;
      if(nfeltval(nf, cbar[1], pr)>0, rspot++; rlis[rspot] = 0);
      cur=1;
      for(k=1, p^(pr[4])-1,
        cur *= resg;
        if(nfeltval(nf, sum(jj=0,poldegree(c), cbar[jj+1]*cur^jj), pr)>0,
          rspot++; rlis[rspot] = cur)),
      /* If p is odd, use another method */
      rrs = reducedroots(nf,pr,c);
      rspot=rrs[1];
      rlis = vector(length(rrs[2]),h,lift(subst(lift(rrs[2][h]),y,x)));
      rlis = vector(length(rlis),h,Mod(rlis[h], g));
      );
    for(k=1, rspot,
      hh= subst(c,x,x*unif+subst(lift(rlis[k]),x,y));
      hh=psharp(hh, p, nf, pr, unif);
      ldeg=reducedpoldegree(hh, nf, pr);
      if(ldeg==1, if(early, return(1), m++));
      if(ldeg>1, cspot++;clis[cspot] = hh));
    );
  m
}

    
/* Counting roots of f an extension of Qp
   unr defines an unramified extension, and is given in terms
   of t.  eis is eisenstein over the unramified extension, so it may
   have t in its coefficients.                                                  */
countrootsdb(p, f, unr1, gg, early=0) =
{       
  local(clis, rlist, cspot, rspot,m, c, eis, hh,ldeg, unif, rrs);
  local(unr);
    
  clis = vector(poldegree(f),h,0);
  rlis = vector(poldegree(f),h,-1);
  cspot=1;
  eis = subst(gg,x,y);  /* in case it was in terms of x */
  if(poldegree(eis,y)==1,eis=y-p);
  unr = subst(unr1,y,t);
  unr = subst(unr,x,t);
  /* My uniformizer is y */

  clis[1]= psharpdb(f, p, eis,unr);
  m=0;
  while(cspot>0,
    c=clis[cspot];
    cspot--;
    /* This is where we find the roots of c modulo p.
    They go into rlis, where the spot is controlled by rspot
    */
    cbar = lift(subst(c,y,0)*Mod(1,p));
    cbar = lift(lift(factorff(cbar, p, unr)))[,1]~; /* irred factors */
    rspot=0;
    for(k=1,length(cbar),
      if(poldegree(cbar[k],x)==1, /* root */
        rspot++;
        rlis[rspot] = lift(lift((-Mod(1,p)*polcoeff(cbar[k],0,x)/polcoeff(cbar[k],1,x))*Mod(1,unr)));
        ));
    /* Now go through the list of roots */
    for(k=1, rspot,
      hh= lift(lift(subst(c,x,x*y+rlis[k])*Mod(1,eis))*Mod(1,unr)); /* ??? */
      hh=psharpdb(hh, p, eis,unr);
      ldeg= poldegree(lift(subst(hh, y,0)*Mod(1,p)),x);
      if(ldeg==1, if(early, return(1), m++));
      if(ldeg>1, cspot++;clis[cspot] = hh));
    );
  m
  }

/* Utility functions for countroots */


/* Take a coefficient which may be in y and t and compute its ord */

upperorddb(ff, p, eram) =
{
  local(ans, k, an2);
  k=0;
  while(polcoeff(ff,k, y)==0, k++);
  ans = coefmaxord(polcoeff(ff,k, y), p)*eram+k;
  while(k<poldegree(ff,y),
    k++;
    if(polcoeff(ff,k, y) != 0,
      an2=coefmaxord(polcoeff(ff,k,y), p)*eram+k;
      ans=min(ans, an2)));
  return(ans);
  }

/* Divide away biggest power of uniformizer */
psharpdb(f, p, eis,unr) =
{
  local(m,j,pc,pe);
  m=999999999999999999999;
  pe=poldegree(eis);
  for(j=0, poldegree(f,x),
    pc=polcoeff(f, j);
    if(pc!=0,  m=min(m, upperorddb(pc, p, pe))));
  if(m>100000, return(f));
  pc=lift(f/Mod(y^m,eis));
  pc=lift(pc*Mod(1,unr));
  return(pc);
  }

/* Take a poly in t and compute the max ord p of all coeffs */

coefmaxord(ff, p) =
{
  local(ans,pc);
  ans=999999999;
  for(j=0,poldegree(ff),
    pc=polcoeff(ff,j,t);
    if(pc != 0, ans=min(ans, valuation(pc,p))));
  return(ans);
  }


psharp(f, p, nf, pr, unif) =
{
  local(m, d, j);

  d = poldegree(f);
  m=999999999999999999999;
  for(j=0, d,
    m=min(m, upperord(polcoeff(f, j), p, nf, pr)));
  if(m>100000, return(f));
  f/unif^m
  }

upperord(ff, p, nf, pr) =
{
  local(j, d, f, g, deg, mycoeff);

  f = lift(ff);
  f = subst(f,y,x);
  g = nf.pol;
  nfeltval(nf, Mod(f, g), pr);
  }


/* Find a generator for the residue field */
nfprimroot(nf, pr) =
{
  local(looking=1, v, g, resord, p, plist, cking, j);
  g=nf.pol;
  p=pr[1];
  if(pr[4]==1, return(lift(znprimroot(p))));
  resord = p^pr[4]-1;
  plist=factor(resord)[,1]~;
  v=vector(poldegree(g),h,0);
  while(looking,
    v=incv(v,p);
    resg=Mod(nf.zk*(v~),g);
    if(nfeltval(nf, resg, pr)==0,
      cking=1; j = 1;
      while(cking && j<=length(plist),
        res2=resg^(resord/plist[j]);
        if(nfeltval(nf, res2-1, pr)>0, cking=0, j++));
      if(cking==1, return(resg))));
} /* Always return from inside */


/* Increment a vector being used as a multi-dimensional counter, each
   component being less than p, and assuming it won't overflow */
incv(v,p) =
{
  local(j=1);
  while(v[j]==p-1, j++);
  v[j] = v[j]+1;
  while(j>1, j--;v[j]=0);
  v
}

/* Degree of this polynomial modulo ideal pr */
reducedpoldegree(pol, nf, pr) =
{
  local(thisdeg);
  thisdeg=poldegree(pol);
  while((thisdeg>-1) && nfeltval(nf,
      Mod(subst(lift(polcoeff(pol, thisdeg)),y,x), nf.pol),
      pr)>0, thisdeg--);
  thisdeg
}