global(MAXPRIME:small=100000);
global(CHECKPRIMEBOUND:small=100);

setMP( maxprime:small ) = {MAXPRIME = maxprime;return;}
setCPB( maxprime:small ) = {CHECKPRIMEBOUND = maxprime;return;}


/* Checks for p giving a degree one prime below check bound.
  Search appropriate for trivial automorphism
*/
check_poly_ngw( poly ) = {
   my(pd:small);
   pd = poldisc(poly):small;

   forprime(p=3,CHECKPRIMEBOUND,
       if( pd % p == 0, next);

       fmp = factormod( poly, p, 1);
       if(fmp[1][1] == 1, return(p));
   );
   return(0);
}

/* finds first degree one prime 
  Search appropriate for trivial automorphism
*/
find_best_poly_ngw( poly ) = {
   my(pd:small,nfd:small,ipr,fmp);
   pd = poldisc(poly):small;
   nfd = nfdisc(poly):small;
   forprime( p=2, MAXPRIME,
       if( nfd % p == 0, next,
           pd % p == 0,
           ipr = idealprimedec(nfinit(poly),p,1);
           if(length(ipr) > 0, return(p);)
           ,
           fmp = factormod( poly, p, 1 );
           if(fmp[1][1] == 1, return(p));
       );
   );
   return(0);
}

/* checks for p giving a non-trivial frobenius of order ord
  NOTE: think about what this means for a given group as it may not give what you want 
        if all order k elements are not conjugate or the field is not galois over Q.
  MOSTLY USEFUL IN C2, S3
*/
check_poly_ng_nt( poly, ord:small ) = {
   my(pd:small,fmp);
   pd = poldisc(poly):small;

   forprime(p=3,CHECKPRIMEBOUND,
       if( pd % p == 0, next);

       fmp = factormod( poly, p, 1);
       if(fmp[length(fmp)][1] == ord, return(p));
   );
   return(0);
}

/* Finds best p giving a non-trivial frobenius of order ord
  NOTE: think about what this means for a given group as it may not give what you want 
        if all order k elements are not conjugate or the field is not galois over Q.
  MOSTLY USEFUL IN C2, S3
*/
find_best_poly_gw_ord( poly, disc:small, ord:small ) = {
   my(pd:small,nfd:small,ipr,fmp);
   pd = poldisc(poly):small;
   nfd = nfdisc(poly):small;
   forprime( p=2, MAXPRIME,
       if( nfd % p == 0, next;
           pd % p == 0,
           ipr = idealprimedec(nfinit(poly),pr,ord);
           for( k:small = 1, length(ipr), 
               if(ipr[k].f == ord, return(p);)
           );
           ,
           fmp = factormod( poly, p, 1 );
           if(fmp[length(fmp)][1] == ord, return(p));
       );
   );
   return(0);
}

/* checks for p giving a a specific frobenius
*/
check_poly_gw( nf, aut ) = {
   my(pd:small,ipr,modP,frobp,autp,twip);
   pd = poldisc(nf.pol):small;

   forprime(p=3,CHECKPRIMEBOUND,
       if( pd % p == 0, next);

       ipr = idealprimedec(nf,p);
       for(j:small=1,length(ipr),
           twip = nfgaloisapply(nf,aut,ipr[j]);
           if( idealval(nf,twip,ipr[j]) != 1, next; );

           modP=nfmodprinit(nf,ipr[j]);
           frobp = nfmodpr(nf,'x,modP)^p;
           autp = nfmodpr(nf,aut,modP);
           if( autp == frobp, return(p));
       );
   );
   return(0);
}

/* searches for p giving a a specific frobenius
  NOTE: it does not search at any place over a prime ramified in L even if some such places are not ramified.
*/
find_best_poly_gw( nf, aut, nfzk ) = {
   my(pd:small,ipr,modP,frobp,autp,twip,IMaut);
   pd = poldisc(nf.pol):small;

   forprime(p=2, MAXPRIME,
       if( nf.disc % p == 0, next);
       ipr = idealprimedec(nf,p);
       if( pd % p == 0,
          for(j:small=1,length(ipr):small,
              twip = nfgaloisapply(nf,aut,ipr[j]);
              if( idealval(nf,twip,ipr[j]) != 1, next; );
              modP=nfmodprinit(nf,ipr[j]);

              my(flag=p);
              for(k:small = 2, length(nfzk):small,
                 IMaut = nfgaloisapply(nf,aut,nfzk[k]);
                 frobp = nfmodpr(nf,nfzk[k],modP)^p;
                 autp = nfmodpr(nf,IMaut,modP);
                 if( autp != frobp, flag=0; break);
              );
              if(flag,return(p));
          );
           ,
          for(j:small=1,length(ipr):small,
              twip = nfgaloisapply(nf,aut,ipr[j]);
              if( idealval(nf,twip,ipr[j]) != 1, next; );
   
              modP=nfmodprinit(nf,ipr[j]);
              frobp = nfmodpr(nf,x,modP)^p;
              autp = nfmodpr(nf,aut,modP);
              if( autp == frobp, return(p));
          );
      );
   );
   return(0);
}