/*-*- compile-command: "
/usr/bin/gcc -c -o check_helper.gp.o -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -I"/usr/local/pari211/include" check_helper.gp.c && /usr/bin/gcc -o check_helper.gp.so -shared -O3 -Wall -fno-strict-aliasing -fomit-frame-pointer -fPIC -Wl,-shared check_helper.gp.o -lc -lm -L/usr/local/pari211/lib -lpari
; -*-*/
#include <pari/pari.h>
/*
GP;install("setMP","vL","setMP","./check_helper.gp.so");
GP;install("check_poly_ngw","lD0,G,L","check_poly_ngw","./check_helper.gp.so");
GP;install("find_best_poly_ngw","lD0,G,L,p","find_best_poly_ngw","./check_helper.gp.so");
GP;install("check_poly_ng_nt","lD0,G,L,L","check_poly_ng_nt","./check_helper.gp.so");
GP;install("find_best_poly_gw_ord","lD0,G,L,L,p","find_best_poly_gw_ord","./check_helper.gp.so");
GP;install("check_poly_gw","lD0,G,D0,G,L","check_poly_gw","./check_helper.gp.so");
GP;install("find_best_poly_gw","lD0,G,D0,G,D0,G,","find_best_poly_gw","./check_helper.gp.so");
*/
void setMP(long maxprime);
long check_poly_ngw(GEN poly, long CHECKPRIMEBOUND);
long find_best_poly_ngw(GEN poly, GEN nfd, long prec);
long find_best_poly_ngw_p(GEN poly, GEN nfd, GEN pfacs, long prec);
long check_poly_ng_nt(GEN poly, long ord, long CHECKPRIMEBOUND);
long find_best_poly_gw_ord(GEN poly, GEN nfd, long ord, long prec);
long find_best_poly_gw_ord_p(GEN poly, GEN nfd, GEN pfacs, long ord, long prec);
long check_poly_gw(GEN nf, GEN aut, long CHECKPRIMEBOUND);
long find_best_poly_gw(GEN nf, GEN aut, GEN nfzk);
/*End of prototype*/

static long MAXPRIME = 1000000;
/*End of global vars*/

void
setMP(long maxprime)	  /* void */
{
  MAXPRIME = maxprime;
}


/* Checks for p giving a degree one prime below check bound.
Search appropriate for trivial automorphism
*/
long
check_poly_ngw(GEN poly, long CHECKPRIMEBOUND)
{
  pari_sp ltop = avma;
  long pd;
  GEN fmp;
  GEN gpd = poldisc0(poly, 0);
  if( is_bigint( gpd ) ) {
    {
      forprime_t iter;	  /* forprime */
      u_forprime_init(&iter, 3, CHECKPRIMEBOUND);;
      {
        pari_sp btop = avma;
        long p;
        while ((p = u_forprime_next(&iter)))
        {
          if (smodis(gpd, p) == 0)
            continue;
          fmp = factormod0(poly, stoi(p), 1);
          if (gequal1(gcoeff(fmp, 1, 1)))
          {
            avma = ltop;
            return p;
          }
          avma = btop;
        }
      }
    }
  } else { 
    pd = itos(gpd);
    {
      forprime_t iter;	  /* forprime */
      u_forprime_init(&iter, 3, CHECKPRIMEBOUND);;
      {
        pari_sp btop = avma;
        long p;
        while ((p = u_forprime_next(&iter)))
        {
          if (smodss(pd, p) == 0)
            continue;
          fmp = factormod0(poly, stoi(p), 1);
          if (gequal1(gcoeff(fmp, 1, 1)))
          {
            avma = ltop;
            return p;
          }
          avma = btop;
        }
      }
    }
  }
  avma = ltop;
  return 0;
}

/* finds first degree one prime 
Search appropriate for trivial automorphism
*/
long
find_best_poly_ngw_p(GEN poly, GEN nfd, GEN pfacs, long prec)
{
  pari_sp ltop = avma;
  long pd;
  GEN ipr, fmp;
  GEN gpd = poldisc0(poly, 0);
  GEN nf = NULL;
  if( is_bigint( gpd ) ) 
  {
    forprime_t iter;	  /* forprime */
    u_forprime_init(&iter, 2, MAXPRIME);;
    {
      pari_sp btop = avma;
      long p;
      while ((p = u_forprime_next(&iter)))
      {
        if (smodis(nfd, p) == 0)
        {
          continue;
        }
        else if( smodis(gpd, p) == 0 )
        {
          if( nf == NULL ) {
              /**/
              if( pfacs != NULL )
              {
                  addprimes(pfacs);

                  GEN tv = cgetg(3, t_VEC);
                  gel(tv,1) = poly;
                  gel(tv,2) = pfacs;
                  nf = nfinit0(tv, 0, prec);

                  if( glength( nfcertify( nf )) > 0 )
                  { 
                     nf = nfinit0(poly, 0, prec);
                  }
                  removeprimes(pfacs);
              }
              else
              {
                 GEN pdf = factorint(gabs(nfd,prec),0);
                 addprimes(gel(pdf,1));

                 GEN tv = cgetg(3, t_VEC);
                 gel(tv,1) = poly;
                 gel(tv,2) = pdf;
                 nf = nfinit0(tv, 0, prec);

                 if( glength( nfcertify( nf )) > 0 )
                 { 
                    nf = nfinit0(poly, 0, prec);
                 }
                 removeprimes(gel(pdf,1));
              }
              nf = gerepileupto(btop, nf );
              btop = avma;
          }
          ipr = idealprimedec_limit_f(nf, stoi(p), 1);
          if (glength(ipr) > 0)
          {
            avma = ltop;
            return p;
          }
        }
        else
        {
          fmp = factormod0(poly, stoi(p), 1);
          if (gequal1(gcoeff(fmp, 1, 1)))
          {
            avma = ltop;
            return p;
          }
        }
        avma = btop;
      }
    }
  } 
  else 
  {
    pd = gtos(gpd);
    forprime_t iter;	  /* forprime */
    u_forprime_init(&iter, 2, MAXPRIME);;
    {
      pari_sp btop = avma;
      long p;
      while ((p = u_forprime_next(&iter)))
      {
        if (smodis(nfd, p) == 0)
        {
          continue;
        }
        else if( smodss(pd, p) == 0 )
        {
          if( nf == NULL ) {
              /**/
              if( pfacs != NULL )
              {
                  addprimes(pfacs);

                  GEN tv = cgetg(3, t_VEC);
                  gel(tv,1) = poly;
                  gel(tv,2) = pfacs;
                  nf = nfinit0(tv, 0, prec);

                  if( glength( nfcertify( nf )) > 0 )
                  { 
                     nf = nfinit0(poly, 0, prec);
                  }
                  removeprimes(pfacs);
              }
              else
              {
                 GEN pdf = factorint(gabs(nfd,prec),0);
                 addprimes(gel(pdf,1));

                 GEN tv = cgetg(3, t_VEC);
                 gel(tv,1) = poly;
                 gel(tv,2) = pdf;
                 nf = nfinit0(tv, 0, prec);

                 if( glength( nfcertify( nf )) > 0 )
                 { 
                    nf = nfinit0(poly, 0, prec);
                 }
                 removeprimes(gel(pdf,1));
              }
              nf = gerepileupto(btop, nf);
              btop = avma;
          }
          ipr = idealprimedec_limit_f(nf, stoi(p), 1);
          if (glength(ipr) > 0)
          {
            avma = ltop;
            return p;
          }
        }
        else
        {
          fmp = factormod0(poly, stoi(p), 1);
          if (gequal1(gcoeff(fmp, 1, 1)))
          {
            avma = ltop;
            return p;
          }
        }
        avma = btop;
      }
    }
  }
  avma = ltop;
  pari_err(e_MISC, "Not enough primes.");
  return 0;
}

long
find_best_poly_ngw(GEN poly, GEN nfd, long prec)
{
   return find_best_poly_ngw_p(poly, nfd, NULL, prec);
}

/* 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
*/
long
check_poly_ng_nt(GEN poly, long ord, long CHECKPRIMEBOUND)
{
  pari_sp ltop = avma;
  GEN fmp;
  GEN gpd = poldisc0(poly, 0);
  {
    forprime_t iter;	  /* forprime */
    u_forprime_init(&iter, 3, CHECKPRIMEBOUND);;
    {
      pari_sp btop = avma;
      long p;
      while ((p = u_forprime_next(&iter)))
      {
        if (smodis(gpd, p) == 0)
          continue;
        fmp = factormod0(poly, stoi(p), 1);
        long nfact = gtos(gel(matsize(fmp), 1));
        if (gequalgs(gcoeff(fmp, nfact, 1), ord))
        {
          avma = ltop;
          return p;
        }
        avma = btop;
      }
    }
  }
  avma = ltop;
  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
*/
long
find_best_poly_gw_ord_p(GEN poly, GEN nfd, GEN pfacs, long ord, long prec)
{
  pari_sp ltop = avma;
  GEN ipr, fmp;
  
  GEN nf = NULL;  
  GEN gpd = poldisc0(poly, 0);
  {
    forprime_t iter;	  /* forprime */
    u_forprime_init(&iter, 2, MAXPRIME);
    {
      pari_sp btop = avma;
      long p;
      while ((p = u_forprime_next(&iter)))
      {
        if (smodis(nfd, p) == 0)
        {
          continue; 
        }
        else if (smodis(gpd, p) == 0)
        {
          long l2;
          if( nf == NULL ) {
              /**/
              if( pfacs != NULL )
              {
                  addprimes(pfacs);

                  GEN tv = cgetg(3, t_VEC);
                  gel(tv,1) = poly;
                  gel(tv,2) = pfacs;
                  nf = nfinit0(tv, 0, prec);

                  if( glength( nfcertify( nf )) > 0 )
                  { 
                     nf = nfinit0(poly, 0, prec);
                  }
                  removeprimes(pfacs);
              }
              else
              {
                 GEN pdf = factorint(gabs(nfd,prec),0);
                 addprimes(gel(pdf,1));

                 GEN tv = cgetg(3, t_VEC);
                 gel(tv,1) = poly;
                 gel(tv,2) = pdf;
                 nf = nfinit0(tv, 0, prec);

                 if( glength( nfcertify( nf )) > 0 )
                 { 
                    nf = nfinit0(poly, 0, prec);
                 }
                 removeprimes(gel(pdf,1));
              }
              

              nf = gerepileupto(btop, nf );
              btop = avma;
          }
          ipr = idealprimedec_limit_f(nf, stoi(p), ord);
          l2 = glength(ipr);
          {
            long k;
            for (k = 1; k <= l2; ++k)
            {
              if (pr_get_f(gel(ipr, k)) == ord)
              {
                avma = ltop;
                return p;
              }
            }
          }
        } 
        else 
        {
          fmp = factormod0(poly, stoi(p), 1);
          long nfact = gtos(gel(matsize(fmp), 1));
          if (gequalgs(gcoeff(fmp, nfact, 1), ord))
          {
            avma = ltop;
            return p;
          }
        }
        avma = btop;
      }
    }
  }
  avma = ltop;
  pari_err(e_MISC, "Not enough primes.");
  return 0;
}


long
find_best_poly_gw_ord(GEN poly, GEN nfd, long ord, long prec)
{
   return find_best_poly_gw_ord_p(poly, nfd, NULL, ord, prec);
}

/* checks for p giving a a specific frobenius
*/
long
check_poly_gw(GEN nf, GEN aut, long CHECKPRIMEBOUND)
{
  pari_sp ltop = avma;
  GEN ipr, modP, frobp, autp, twip;
  GEN x = pol_x(0); /* TODO global ?? */
  GEN gpd = poldisc0(member_pol(nf), 0);
  {
    forprime_t iter;	  /* forprime */
    u_forprime_init(&iter, 3, CHECKPRIMEBOUND);;
    {
      pari_sp btop = avma;
      long p;
      while ((p = u_forprime_next(&iter)))
      {
        {
          long l1;
          if (smodis(gpd, p) == 0)
            continue;
          ipr = idealprimedec(nf, stoi(p));
          l1 = glength(ipr);
          {
            long j;
            for (j = 1; j <= l1; ++j)
            {
              twip = galoisapply(nf, aut, gel(ipr, j));
              if (!gequalgs(gpidealval(nf, twip, gel(ipr, j)), 1))
                continue;
              modP = nfmodprinit(nf, gel(ipr, j));
              frobp = gpowgs(nfmodpr(nf, x, modP), p);
              autp = nfmodpr(nf, aut, modP);
              if (gequal(autp, frobp))
              {
                avma = ltop;
                return p;
              }
            }
          }
        }
        avma = btop;
      }
    }
  }
  avma = ltop;
  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.
*/
long
find_best_poly_gw(GEN nf, GEN aut, GEN nfzk)
{
  pari_sp ltop = avma;
  GEN ipr, modP, frobp, autp, twip, IMaut;
  GEN x = pol_x(0);
  GEN gpd = poldisc0(member_pol(nf), 0);
  {
    forprime_t iter;	  /* forprime */
    u_forprime_init(&iter, 2, MAXPRIME);;
    {
      pari_sp btop = avma;
      long p;
      while ((p = u_forprime_next(&iter)))
      {
        if (gequal0(gmodgs(member_disc(nf), p)))
          continue;
        ipr = idealprimedec(nf, stoi(p));
        if (smodis(gpd, p) == 0) 
        {
          long l1;
          l1 = glength(ipr);
          {
            long j;
            for (j = 1; j <= l1; ++j)
            {
              {
                long l2;
                twip = galoisapply(nf, aut, gel(ipr, j));
                if (!gequalgs(gpidealval(nf, twip, gel(ipr, j)), 1))
                  continue;
                modP = nfmodprinit(nf, gel(ipr, j));
                long flag = 1;
                l2 = glength(nfzk);
                {
                  long k;
                  for (k = 2; k <= l2; ++k)
                  {
                    IMaut = galoisapply(nf, aut, gel(nfzk, k));
                    frobp = gpowgs(nfmodpr(nf, gel(nfzk, k), modP), p);
                    autp = nfmodpr(nf, IMaut, modP);
                    if (!gequal(autp, frobp))
                    {
                      flag = 0;
                      break;
                    }
                  }
                }
                if (flag == 1)
                {
                  avma = ltop;
                  return p;
                }
              }
            }
          }
        }
        else
        {
          long l4;
          l4 = glength(ipr);
          {
            long j;
            for (j = 1; j <= l4; ++j)
            {
              twip = galoisapply(nf, aut, gel(ipr, j));
              if (!gequalgs(gpidealval(nf, twip, gel(ipr, j)), 1))
                continue;
              modP = nfmodprinit(nf, gel(ipr, j));
              frobp = gpowgs(nfmodpr(nf, x, modP), p);
              autp = nfmodpr(nf, aut, modP);
              if (gequal(autp, frobp))
              {
                avma = ltop;
                return p;
              }
            }
          }
        }
        avma = btop;
      }
    }
  }
  avma = ltop;
  pari_err(e_MISC, "Not enough primes.");
  return 0;
}


void
find_best_file_ng(char *fname, long prec)         /* void */
{
  pari_sp ltop = avma;
  GEN nextpoly;
  long fn = gp_fileopen(fname, "r");
  {
    pari_sp btop = avma;
    while (!gequal0(nextpoly = geval(gp_fileread(fn))))
    {
      if( glength( nextpoly ) > 2 )
      {
         pari_printf("[%Ps,%Ps,%ld]\n", gel(nextpoly, 2), gel(nextpoly, 1), find_best_poly_ngw_p(gel(nextpoly, 2), gel(nextpoly, 1), gel(nextpoly, 3), prec));
      }
      else
      {
         pari_printf("[%Ps,%Ps,%ld]\n", gel(nextpoly, 2), gel(nextpoly, 1), find_best_poly_ngw(gel(nextpoly, 2), gel(nextpoly, 1), prec));
      }
      avma = btop;
    }
  }
  gp_fileclose(fn);
  avma = ltop;
  return;
}


void
find_best_file_ng_of(char *fname, long ofn, long prec)    /* void */
{
  pari_sp ltop = avma;
  GEN nextpoly;
  long fn = gp_fileopen(fname, "r");
  {
    pari_sp btop = avma;
    while (!gequal0(nextpoly = geval(gp_fileread(fn))))
    {
      long p;
      if( glength( nextpoly ) > 2 )
      {
         p = find_best_poly_ngw_p(gel(nextpoly, 2), gel(nextpoly, 1), gel(nextpoly, 3), prec);
      }
      else
      {
         p = find_best_poly_ngw(gel(nextpoly, 2), gel(nextpoly, 1), prec);
      }
      char * str = stack_sprintf("[%Ps,%Ps,%ld]",gel(nextpoly, 2),gel(nextpoly, 1),p);
      gp_filewrite(ofn, str);
      avma = btop;
    }
  }
  gp_fileclose(fn);
  avma = ltop;
  return;
}

void
find_best_file_allaut(char *fname, long num_auts, long prec)         /* void */
{
  pari_sp ltop = avma;
  long fn;
  GEN nf, nfzk, auts, nextpoly;
  fn = gp_fileopen(fname, "r");
  {
    pari_sp btop = avma;
    while (!gequal0(nextpoly = geval(gp_fileread(fn))))
    {
      {
        long l2;
        /* */
        if( glength( nextpoly ) > 2 )
        {
           addprimes(gel(nextpoly, 3));
           GEN p1 = cgetg(3, t_VEC);
           gel(p1, 1) = gcopy(gel(nextpoly, 2));
           gel(p1, 2) = gel(nextpoly, 3);;
           nf = nfinit0(p1, 0, prec);

           if( glength( nfcertify( nf )) > 0 )
           { 
              nf = nfinit0(gel(nextpoly, 2), 0, prec);
           }
           removeprimes(gel(nextpoly, 3));
        }
        else
        {
           GEN pdf = factorint(gabs(gel(nextpoly, 1),prec),0);
           addprimes(gel(pdf,1));

           GEN p1 = cgetg(3, t_VEC);
           gel(p1, 1) = gcopy(gel(nextpoly, 2));
           gel(p1, 2) = pdf;
           nf = nfinit0(p1, 0, prec);

           if( glength( nfcertify( nf )) > 0 )
           { 
              nf = nfinit0(gel(nextpoly, 2), 0, prec);
           }
           removeprimes(gel(pdf,1));
        }
        auts = galoisconj0(nf, 0, NULL, prec);
        if (glength(auts) != num_auts)
          pari_err(e_MISC, "Wrong Number of Auts");
        nfzk = member_zk(nf);
        l2 = glength(auts);
        {
          long i;
          for (i = 1; i <= l2; ++i)
          {
            pari_printf("[%Ps,%Ps,%ld,%Ps]\n", gel(nextpoly, 2), gel(nextpoly, 1), find_best_poly_gw(nf, gel(auts, i), nfzk), gel(auts, i));
          }
        }
      }
      avma=btop;
    }
  }
  gp_fileclose(fn);
  avma = ltop;
  return;
}

void
find_best_file_allaut_of(char * fname, long num_auts, long ofn, long prec)     /* void */
{
  pari_sp ltop = avma;
  long fn;
  GEN nf, nfzk, auts, nextpoly;
  fn = gp_fileopen(fname, "r");
  {
    pari_sp btop = avma;
    while (!gequal0(nextpoly = geval(gp_fileread(fn))))
    {
      {
        long l2;
       /* 
      */
        if( glength( nextpoly ) > 2 )
        {
           //pari_printf("%Ps\n", nextpoly );
           addprimes(gel(nextpoly, 3));

           GEN p1 = cgetg(3, t_VEC);
           gel(p1, 1) = gel(nextpoly, 2);
           gel(p1, 2) = gel(nextpoly, 3);
           nf = nfinit0(p1, 0, prec);
   
           if( glength( nfcertify( nf )) > 0 )
           { 
              nf = nfinit0(gel(nextpoly, 2), 0, prec);
           }
           removeprimes(gel(nextpoly, 3));
        }
        else
        {
           GEN pdf = factorint(gabs(gel(nextpoly, 1),prec),0);
           addprimes(gel(pdf,1));
           GEN p1 = cgetg(3, t_VEC);
           gel(p1, 1) = gel(nextpoly, 2);
           gel(p1, 2) = pdf;
           nf = nfinit0(p1, 0, prec);
   
           if( glength( nfcertify( nf )) > 0 )
           { 
              nf = nfinit0(gel(nextpoly, 2), 0, prec);
           }
           removeprimes(gel(pdf,1));
        }
        auts = galoisconj0(nf, 0, NULL, prec);
        if (glength(auts) != num_auts)
        { 
          pari_printf("%Ps %ld %ld\n", nextpoly, num_auts, glength(auts));
          pari_err(e_MISC, "Wrong Number of Auts %Ps", nextpoly);
        }
        nfzk = member_zk(nf);
        l2 = glength(auts);
        {
          long i;
          for (i = 1; i <= l2; ++i)
          {
            long p = find_best_poly_gw(nf, gel(auts, i), nfzk);
            char * str = stack_sprintf("[%Ps,%Ps,%ld,%Ps]",gel(nextpoly, 2),gel(nextpoly, 1),p, gel(auts, i));
            gp_filewrite(ofn, str);
          }
        }
      }
      avma = btop;
    }
  }
  gp_fileclose(fn);
  avma = ltop;
  return;
}

void
find_best_file_C2(char * fname, long prec)   /* void */
{
  pari_sp ltop = avma;
  GEN nextpoly;
  long fn = gp_fileopen(fname, "r");
  {
    pari_sp btop = avma;
    while (!gequal0(nextpoly = geval(gp_fileread(fn))))
    {
      pari_printf("[%Ps,%Ps,%ld,%ld]\n", gel(nextpoly, 2), gel(nextpoly, 1), find_best_poly_ngw(gel(nextpoly, 2), gel(nextpoly, 1), prec), find_best_poly_gw_ord(gel(nextpoly, 2), gel(nextpoly, 1), 2, prec));
      avma = btop;
    }
  }
  gp_fileclose(fn);
  avma = ltop;
  return;
}

void
find_best_file_C2_of(char * fname, long ofn, long prec)       /* void */
{
  pari_sp ltop = avma;
  long fn;
  GEN nextpoly = gen_0;
  fn = gp_fileopen(fname, "r");
  {
    pari_sp btop = avma;
    while (!gequal0(nextpoly = geval(gp_fileread(fn))))
    {
      long p1 = find_best_poly_ngw(gel(nextpoly, 2), gel(nextpoly, 1), prec);
      long p2 = find_best_poly_gw_ord(gel(nextpoly, 2), gel(nextpoly, 1), 2, prec);
      char * str = stack_sprintf("[%Ps,%Ps,%ld,%ld]",gel(nextpoly, 2),gel(nextpoly, 1),p1,p2);
      gp_filewrite(ofn, str);
      avma = btop;
    }
  }
  gp_fileclose(fn);
  avma = ltop;
  return;
}


void
iter_best_C2_of(long end, long ofn, long prec)
{
  pari_sp ltop = avma;
  GEN xx = pol_x(0);
  {
    pari_sp btop = avma;
    long i;
    for (i = 1; i <= end; ++i)
    {
      if (issquarefree(stoi(1 - (4*i)))) {
        GEN pol = gaddgs(gsub(gsqr(xx), xx), i);
        GEN dis = stoi(1 - (4*i));
        long p1 = find_best_poly_ngw(pol, dis, prec);
        long p2 = find_best_poly_gw_ord(pol, dis, 2, prec);
        char * str = stack_sprintf("[%Ps,%Ps,%ld,%ld]",pol,dis,p1,p2);
        gp_filewrite(ofn, str);
      }
      if (issquarefree(stoi(1 + (4*i)))){
        GEN pol = gsubgs(gsub(gsqr(xx), xx), i);
        GEN dis = stoi(1 + (4*i));
        long p1 = find_best_poly_ngw(pol, dis, prec);
        long p2 = find_best_poly_gw_ord(pol, dis, 2, prec);
        char * str = stack_sprintf("[%Ps,%Ps,%ld,%ld]",pol,dis,p1,p2);
        gp_filewrite(ofn, str);
      }
      if (issquarefree(stoi(i)))
      {
        if (smodss(i, 4) == 3){
          GEN pol = gsubgs(gsqr(xx), i);
          GEN dis = stoi(4*i);
          long p1 = find_best_poly_ngw(pol, dis, prec);
          long p2 = find_best_poly_gw_ord(pol, dis, 2, prec);
          char * str = stack_sprintf("[%Ps,%Ps,%ld,%ld]",pol,dis,p1,p2);
          gp_filewrite(ofn, str);
        }
        if (smodss(i, 4) == 1){
          GEN pol = gaddgs(gsqr(xx), i);
          GEN dis = stoi(-4*i);
          long p1 = find_best_poly_ngw(pol, dis, prec);
          long p2 = find_best_poly_gw_ord(pol, dis, 2, prec);
          char * str = stack_sprintf("[%Ps,%Ps,%ld,%ld]",pol,dis,p1,p2);
          gp_filewrite(ofn, str);
        }
        if (smodss(i, 4) == 2){
          GEN pol = gaddgs(gsqr(xx), i);
          GEN dis = stoi(-4*i);
          long p1 = find_best_poly_ngw(pol, dis, prec);
          long p2 = find_best_poly_gw_ord(pol, dis, 2, prec);
          char * str = stack_sprintf("[%Ps,%Ps,%ld,%ld]",pol,dis,p1,p2);
          gp_filewrite(ofn, str);
        }
        if (smodss(i, 4) == 2){
          GEN pol = gsubgs(gsqr(xx), i);
          GEN dis = stoi(4*i);
          long p1 = find_best_poly_ngw(pol, dis, prec);
          long p2 = find_best_poly_gw_ord(pol, dis, 2, prec);
          char * str = stack_sprintf("[%Ps,%Ps,%ld,%ld]",pol,dis,p1,p2);
          gp_filewrite(ofn, str);
        }
      }
      avma = btop;
    }
  }
  avma = ltop;
  return;
}