[BACK]Return to base1.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/base1.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:01 version 1.2, 2002/09/11 07:26:48
Line 20  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
Line 20  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /**************************************************************/  /**************************************************************/
 #include "pari.h"  #include "pari.h"
 #include "parinf.h"  #include "parinf.h"
 GEN idealhermite_aux(GEN nf, GEN x);  
   
   int new_galois_format = 0;
   
   extern GEN R_from_QR(GEN x, long prec);
   extern GEN to_polmod(GEN x, GEN mod);
   extern GEN roots_to_pol_r1r2(GEN a, long r1, long v);
   extern GEN idealhermite_aux(GEN nf, GEN x);
   extern GEN cauchy_bound(GEN p);
   extern GEN galoisbig(GEN x, long prec);
   extern GEN lllfp_marked(int M, GEN x, long D, long flag, long prec, int gram);
   extern GEN lllint_marked(long M, GEN x, long D, int g, GEN *h, GEN *fl, GEN *B);
   extern GEN mulmat_pol(GEN A, GEN x);
   extern ulong smulss(ulong x, ulong y, ulong *rem);
   
 void  void
 checkrnf(GEN rnf)  checkrnf(GEN rnf)
 {  {
   if (typ(rnf)!=t_VEC) err(idealer1);    if (typ(rnf)!=t_VEC || lg(rnf)!=12) err(idealer1);
   if (lg(rnf)!=12) err(idealer1);  
 }  }
   
 GEN  GEN
 checkbnf(GEN bnf)  _checkbnf(GEN bnf)
 {  {
   if (typ(bnf)!=t_VEC) err(idealer1);    if (typ(bnf) == t_VEC)
   switch (lg(bnf))      switch (lg(bnf))
       {
         case 11: return bnf;
         case 7:  return checkbnf((GEN)bnf[1]);
   
         case 3:
           if (typ(bnf[2])==t_POLMOD)
             return checkbnf((GEN)bnf[1]);
       }
     return NULL;
   }
   
   GEN
   _checknf(GEN nf)
   {
     if (typ(nf)==t_VEC)
       switch(lg(nf))
       {
         case 10: return nf;
         case 11: return checknf((GEN)nf[7]);
         case 7:  return checknf((GEN)nf[1]);
         case 3: if (typ(nf[2]) == t_POLMOD) return checknf((GEN)nf[1]);
       }
     return NULL;
   }
   
   GEN
   checkbnf(GEN x)
   {
     GEN bnf = _checkbnf(x);
     if (!bnf)
   {    {
     case 11: return bnf;      if (_checknf(x)) err(talker,"please apply bnfinit first");
     case 10:      err(idealer1);
       if (typ(bnf[1])==t_POL)  
         err(talker,"please apply bnfinit first");  
       break;  
     case 7:  
       return checkbnf((GEN)bnf[1]);  
   
     case 3:  
       if (typ(bnf[2])==t_POLMOD)  
         return checkbnf((GEN)bnf[1]);  
   }    }
   err(idealer1);    return bnf;
   return NULL; /* not reached */  
 }  }
   
 GEN  GEN
Line 61  checkbnf_discard(GEN bnf)
Line 92  checkbnf_discard(GEN bnf)
 }  }
   
 GEN  GEN
 checknf(GEN nf)  checknf(GEN x)
 {  {
   if (typ(nf)==t_POL) err(talker,"please apply nfinit first");    GEN nf = _checknf(x);
   if (typ(nf)!=t_VEC) err(idealer1);    if (!nf)
   switch(lg(nf))  
   {    {
     case 10: return nf;      if (typ(x)==t_POL) err(talker,"please apply nfinit first");
     case 11: return checknf((GEN)nf[7]);      err(idealer1);
     case 7:  return checknf((GEN)nf[1]);  
     case 3: if (typ(nf[2]) == t_POLMOD) return checknf((GEN)nf[1]);  
   }    }
   err(idealer1);    return nf;
   return NULL; /* not reached */  
 }  }
   
 void  void
Line 126  checkprimeid(GEN id)
Line 153  checkprimeid(GEN id)
 }  }
   
 void  void
 checkprhall(GEN prhall)  
 {  
   if (typ(prhall) != t_VEC || lg(prhall) != 3 || typ(prhall[1]) != t_MAT)  
     err(talker,"incorrect prhall format");  
 }  
   
 void  
 check_pol_int(GEN x, char *s)  check_pol_int(GEN x, char *s)
 {  {
   long k = lgef(x)-1;    long k = lgef(x)-1;
Line 248  transroot(GEN x, int i, int j)
Line 268  transroot(GEN x, int i, int j)
 GEN  GEN
 tschirnhaus(GEN x)  tschirnhaus(GEN x)
 {  {
   long a, v = varn(x), av = avma;    gpmem_t av = avma, av2;
     long a, v = varn(x);
   GEN u, p1 = cgetg(5,t_POL);    GEN u, p1 = cgetg(5,t_POL);
   
   if (typ(x)!=t_POL) err(notpoler,"tschirnhaus");    if (typ(x)!=t_POL) err(notpoler,"tschirnhaus");
Line 260  tschirnhaus(GEN x)
Line 281  tschirnhaus(GEN x)
     a = mymyrand() >> randshift; if (a==0) a =1; p1[4]=lstoi(a);      a = mymyrand() >> randshift; if (a==0) a =1; p1[4]=lstoi(a);
     a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[3]=lstoi(a);      a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[3]=lstoi(a);
     a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[2]=lstoi(a);      a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[2]=lstoi(a);
     u = caract2(x,p1,v); a=avma;      u = caract2(x,p1,v); av2=avma;
   }    }
   while (lgef(srgcd(u,derivpol(u))) > 3); /* while u not separable */    while (lgef(srgcd(u,derivpol(u))) > 3); /* while u not separable */
   if (DEBUGLEVEL>1)    if (DEBUGLEVEL>1)
     fprintferr("Tschirnhaus transform. New pol: %Z",u);      fprintferr("Tschirnhaus transform. New pol: %Z",u);
   avma=a; return gerepileupto(av,u);    avma=av2; return gerepileupto(av,u);
 }  }
 #undef randshift  #undef randshift
   
Line 284  gpolcomp(GEN p1, GEN p2)
Line 305  gpolcomp(GEN p1, GEN p2)
   return 0;    return 0;
 }  }
   
 /* assume pol in Z[X] */  /* assume pol in Z[X]. Find C, L in Z such that POL = C pol(x / L) monic in Z[X].
    * Return POL and set *ptlead = L */
 GEN  GEN
 primitive_pol_to_monic(GEN pol, GEN *ptlead)  primitive_pol_to_monic(GEN pol, GEN *ptlead)
 {  {
Line 342  get_F4(GEN x)
Line 364  get_F4(GEN x)
   return p1;    return p1;
 }  }
   
 GEN galoisbig(GEN x, long prec);  static GEN
 GEN roots_to_pol(GEN a, long v);  roots_to_ZX(GEN z, long *e)
   {
     GEN a = roots_to_pol(z,0);
     GEN b = grndtoi(greal(a),e);
     long e1 = gexpo(gimag(a));
     if (e1 > *e) *e = e1;
     return b;
   }
   
   static GEN
   _res(long n, long s, long k)
   {
     GEN y = cgetg(4, t_VEC);
     y[1] = lstoi(n);
     y[2] = lstoi(s);
     if (!new_galois_format) k = (n == 6 && (k == 6 || k == 2))? 2: 1;
     y[3] = lstoi(k); return y;
   }
   
 GEN  GEN
 galois(GEN x, long prec)  galois(GEN x, long prec)
 {  {
   long av=avma,av1,i,j,k,n,f,l,l2,e,e1,pr,ind;    gpmem_t av = avma, av1;
   GEN x1,p1,p2,p3,p4,p5,p6,w,y,z,ee;    long i,j,k,n,f,l,l2,e,e1,pr,ind;
     GEN x1,p1,p2,p3,p4,p5,w,z,ee;
   static int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3};    static int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3};
   static int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4,    static int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4,
                        1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5,                         1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5,
Line 357  galois(GEN x, long prec)
Line 397  galois(GEN x, long prec)
   if (typ(x)!=t_POL) err(notpoler,"galois");    if (typ(x)!=t_POL) err(notpoler,"galois");
   n=degpol(x); if (n<=0) err(constpoler,"galois");    n=degpol(x); if (n<=0) err(constpoler,"galois");
   if (n>11) err(impl,"galois of degree higher than 11");    if (n>11) err(impl,"galois of degree higher than 11");
   x = gdiv(x,content(x));    x = primpart(x);
   check_pol_int(x, "galois");    check_pol_int(x, "galois");
   if (gisirreducible(x) != gun)    if (gisirreducible(x) != gun)
     err(impl,"galois of reducible polynomial");      err(impl,"galois of reducible polynomial");
   
   if (n<4)    if (n<4)
   {    {
     if (n<3)      if (n == 1) { avma = av; return _res(1, 1,1); }
     {      if (n == 2) { avma = av; return _res(2,-1,1); }
       avma=av; y=cgetg(4,t_VEC);      /* n = 3 */
       y[1] = (n==1)? un: deux;      f = carreparfait(ZX_disc(x));
       y[2]=lnegi(gun);      avma = av;
     }      return f? _res(3,1,1): _res(6,-1,2);
     else /* n=3 */  
     {  
       f=carreparfait(discsr(x));  
       avma=av; y=cgetg(4,t_VEC);  
       if (f) { y[1]=lstoi(3); y[2]=un; }  
       else   { y[1]=lstoi(6); y[2]=lnegi(gun); }  
     }  
     y[3]=un; return y;  
   }    }
   x1 = x = primitive_pol_to_monic(x,NULL); av1=avma;    x1 = x = primitive_pol_to_monic(x,NULL); av1=avma;
   if (n>7) return galoisbig(x,prec);    if (n > 7) return galoisbig(x,prec);
   for(;;)    for(;;)
   {    {
       GEN cb = cauchy_bound(x);
     switch(n)      switch(n)
     {      {
       case 4: z = cgetg(7,t_VEC);        case 4: z = cgetg(7,t_VEC);
           prec = DEFAULTPREC + (long)((gexpo(cb)*18. / BITS_IN_LONG));
         for(;;)          for(;;)
         {          {
           p1=roots(x,prec);            p1=roots(x,prec);
Line 395  galois(GEN x, long prec)
Line 429  galois(GEN x, long prec)
           z[4] = (long)get_F4(transroot(p1,1,4));            z[4] = (long)get_F4(transroot(p1,1,4));
           z[5] = (long)get_F4(transroot(p1,2,3));            z[5] = (long)get_F4(transroot(p1,2,3));
           z[6] = (long)get_F4(transroot(p1,3,4));            z[6] = (long)get_F4(transroot(p1,3,4));
           p4 = roots_to_pol(z,0);            p5 = roots_to_ZX(z, &e);
           p5=grndtoi(greal(p4),&e);  
           e1=gexpo(gimag(p4)); if (e1>e) e=e1;  
           if (e <= -10) break;            if (e <= -10) break;
           prec = (prec<<1)-2;            prec = (prec<<1)-2;
         }          }
         p6 = modulargcd(p5, derivpol(p5));          if (!ZX_is_squarefree(p5)) goto tchi;
         if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;  
         p2 = (GEN)factor(p5)[1];          p2 = (GEN)factor(p5)[1];
         switch(lg(p2)-1)          switch(lg(p2)-1)
         {          {
           case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);            case 1: f = carreparfait(ZX_disc(x)); avma = av;
             y[3]=un;              return f? _res(12,1,4): _res(24,-1,5);
             if (f) { y[2]=un; y[1]=lstoi(12); return y; }  
             y[2]=lnegi(gun); y[1]=lstoi(24); return y;  
   
           case 2: avma=av; y=cgetg(4,t_VEC);            case 2: avma = av; return _res(8,-1,3);
             y[3]=un; y[2]=lnegi(gun); y[1]=lstoi(8); return y;  
   
           case 3: avma=av; y=cgetg(4,t_VEC);            case 3: avma = av;
             y[1]=lstoi(4); y[3]=un;              return (degpol(p2[1])==2)? _res(4,1,2): _res(4,-1,1);
             y[2] = (lgef(p2[1])==5)? un: lnegi(gun);  
             return y;  
   
           default: err(bugparier,"galois (bug1)");            default: err(bugparier,"galois (bug1)");
         }          }
   
       case 5: z = cgetg(7,t_VEC);        case 5: z = cgetg(7,t_VEC);
         ee = new_chunk(7); w = new_chunk(7);          ee= cgetg(7,t_VECSMALL);
           w = cgetg(7,t_VECSMALL);
           prec = DEFAULTPREC + (long)((gexpo(cb)*21. / BITS_IN_LONG));
         for(;;)          for(;;)
         {          {
           for(;;)            for(;;)
Line 443  galois(GEN x, long prec)
Line 471  galois(GEN x, long prec)
               e1 = gexpo(gimag(p3)); if (e1>e) e=e1;                e1 = gexpo(gimag(p3)); if (e1>e) e=e1;
               ee[l]=e; z[l] = (long)p3;                ee[l]=e; z[l] = (long)p3;
             }              }
             p4 = roots_to_pol(z,0);              p5 = roots_to_ZX(z, &e);
             p5=grndtoi(greal(p4),&e);  
             e1 = gexpo(gimag(p4)); if (e1>e) e=e1;  
             if (e <= -10) break;              if (e <= -10) break;
             prec = (prec<<1)-2;              prec = (prec<<1)-2;
           }            }
           p6 = modulargcd(p5,derivpol(p5));            if (!ZX_is_squarefree(p5)) goto tchi;
           if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;  
           p3=(GEN)factor(p5)[1];            p3=(GEN)factor(p5)[1];
           f=carreparfait(discsr(x));            f=carreparfait(ZX_disc(x));
           if (lg(p3)-1==1)            if (lg(p3)-1==1)
           {            {
             avma=av; y=cgetg(4,t_VEC); y[3]=un;              avma = av;
             if (f) { y[2]=un; y[1]=lstoi(60); return y; }              return f? _res(60,1,4): _res(120,-1,5);
             else { y[2]=lneg(gun); y[1]=lstoi(120); return y; }  
           }            }
           if (!f)            if (!f) { avma = av; return _res(20,-1,3); }
           {  
             avma=av; y=cgetg(4,t_VEC);  
             y[3]=un; y[2]=lneg(gun); y[1]=lstoi(20); return y;  
           }  
           pr = - (bit_accuracy(prec) >> 1);            pr = - (bit_accuracy(prec) >> 1);
           for (l=1; l<=6; l++)            for (l=1; l<=6; l++)
             if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break;              if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break;
Line 472  galois(GEN x, long prec)
Line 493  galois(GEN x, long prec)
           p3=gzero;            p3=gzero;
           for (i=1; i<=5; i++)            for (i=1; i<=5; i++)
           {            {
             j=(i%5)+1;              j = i+1; if (j>5) j -= 5;
             p3=gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]),              p3 = gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]),
                             gsub((GEN)p2[j],(GEN)p2[i])));                                gsub((GEN)p2[j],(GEN)p2[i])));
           }            }
           p5=gsqr(p3); p4=grndtoi(greal(p5),&e);            p5=gsqr(p3); p4=grndtoi(greal(p5),&e);
           e1 = gexpo(gimag(p5)); if (e1>e) e=e1;            e1 = gexpo(gimag(p5)); if (e1>e) e=e1;
           if (e <= -10)            if (e <= -10)
           {            {
             if (gcmp0(p4)) goto tchi;              if (gcmp0(p4)) goto tchi;
             f=carreparfait(p4); avma=av; y=cgetg(4,t_VEC);              f = carreparfait(p4); avma = av;
             y[3]=y[2]=un; y[1]=lstoi(f?5:10);              return f? _res(5,1,1): _res(10,1,2);
             return y;  
           }            }
           prec=(prec<<1)-2;            prec=(prec<<1)-2;
         }          }
   
       case 6: z = cgetg(7, t_VEC);        case 6: z = cgetg(7, t_VEC);
           prec = DEFAULTPREC + (long)((gexpo(cb)*42. / BITS_IN_LONG));
         for(;;)          for(;;)
         {          {
           for(;;)            for(;;)
Line 502  galois(GEN x, long prec)
Line 523  galois(GEN x, long prec)
               {                {
                 p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]),                  p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]),
                         gmul((GEN)p2[ind6[k+2]],(GEN)p2[ind6[k+3]]));                          gmul((GEN)p2[ind6[k+2]],(GEN)p2[ind6[k+3]]));
                 p3=gadd(p3,gmul(gsqr(gmul((GEN)p2[i],(GEN)p2[j])),p5)); k+=4;                  p3=gadd(p3,gmul(gsqr(gmul((GEN)p2[i],(GEN)p2[j])),p5));
                   k += 4;
               }                }
               z[l] = (long)p3;                z[l] = (long)p3;
             }              }
             p4=roots_to_pol(z,0);              p5 = roots_to_ZX(z, &e);
             p5=grndtoi(greal(p4),&e);  
             e1 = gexpo(gimag(p4)); if (e1>e) e=e1;  
             if (e <= -10) break;              if (e <= -10) break;
             prec=(prec<<1)-2;              prec=(prec<<1)-2;
           }            }
           p6 = modulargcd(p5,derivpol(p5));            if (!ZX_is_squarefree(p5)) goto tchi;
           if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;  
           p2=(GEN)factor(p5)[1];            p2=(GEN)factor(p5)[1];
           switch(lg(p2)-1)            switch(lg(p2)-1)
           {            {
Line 530  galois(GEN x, long prec)
Line 549  galois(GEN x, long prec)
                           gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6]));                            gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6]));
                   z[++ind] = (long)p3;                    z[++ind] = (long)p3;
                 }                  }
               p4 = roots_to_pol(z,0);                p5 = roots_to_ZX(z, &e);
               p5=grndtoi(greal(p4),&e);  
               e1 = gexpo(gimag(p4)); if (e1>e) e=e1;  
               if (e <= -10)                if (e <= -10)
               {                {
                 p6 = modulargcd(p5,derivpol(p5));                  if (!ZX_is_squarefree(p5)) goto tchi;
                 if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;                  p2 = (GEN)factor(p5)[1];
                 p2=(GEN)factor(p5)[1];                  f = carreparfait(ZX_disc(x));
                 f=carreparfait(discsr(x));                  avma = av;
                 avma=av; y=cgetg(4,t_VEC); y[3]=un;  
                 if (lg(p2)-1==1)                  if (lg(p2)-1==1)
                 {                    return f? _res(360,1,15): _res(720,-1,16);
                   if (f) { y[2]=un; y[1]=lstoi(360); }  
                   else { y[2]=lnegi(gun); y[1]=lstoi(720); }  
                 }  
                 else                  else
                 {                    return f? _res(36,1,10): _res(72,-1,13);
                   if (f) { y[2]=un; y[1]=lstoi(36); }  
                   else { y[2]=lnegi(gun); y[1]=lstoi(72); }  
                 }  
                 return y;  
               }                }
               prec=(prec<<1)-2; break;                prec=(prec<<1)-2; break;
   
             case 2: l2=degpol(p2[1]); if (l2>3) l2=6-l2;              case 2: l2=degpol(p2[1]); if (l2>3) l2=6-l2;
               switch(l2)                switch(l2)
               {                {
                 case 1: f=carreparfait(discsr(x));                  case 1: f = carreparfait(ZX_disc(x));
                   avma=av; y=cgetg(4,t_VEC); y[3]=un;                    avma = av;
                   if (f) { y[2]=un; y[1]=lstoi(60); }                    return f? _res(60,1,12): _res(120,-1,14);
                   else { y[2]=lneg(gun); y[1]=lstoi(120); }                  case 2: f = carreparfait(ZX_disc(x));
                   return y;                    if (f) { avma = av; return _res(24,1,7); }
                 case 2: f=carreparfait(discsr(x));                    p3 = (degpol(p2[1])==2)? (GEN)p2[2]: (GEN)p2[1];
                   if (f)                    f = carreparfait(ZX_disc(p3));
                   {                    avma = av;
                     avma=av; y=cgetg(4,t_VEC);                    return f? _res(24,-1,6): _res(48,-1,11);
                     y[3]=y[2]=un; y[1]=lstoi(24);                  case 3: f = carreparfait(ZX_disc((GEN)p2[1]))
                   }                           || carreparfait(ZX_disc((GEN)p2[2]));
                   else                    avma = av;
                   {                    return f? _res(18,-1,5): _res(36,-1,9);
                     p3=(lgef(p2[1])==5) ? (GEN)p2[2]:(GEN)p2[1];  
                     f=carreparfait(discsr(p3));  
                     avma=av; y=cgetg(4,t_VEC); y[2]=lneg(gun);  
                     if (f) { y[1]=lstoi(24); y[3]=deux; }  
                     else { y[1]=lstoi(48); y[3]=un; }  
                   }  
                   return y;  
                 case 3: f=carreparfait(discsr((GEN)p2[1]))  
                        || carreparfait(discsr((GEN)p2[2]));  
                   avma=av; y=cgetg(4,t_VEC);  
                   y[3]=un; y[2]=lneg(gun); y[1]=lstoi(f? 18: 36);  
                   return y;  
               }                }
             case 3:              case 3:
               for (l2=1; l2<=3; l2++)                for (l2=1; l2<=3; l2++)
                 if (lgef(p2[l2])>=6) p3=(GEN)p2[l2];                  if (degpol(p2[l2]) >= 3) p3 = (GEN)p2[l2];
               if (lgef(p3)==6)                if (degpol(p3) == 3)
               {                {
                 f=carreparfait(discsr(p3)); avma=av; y=cgetg(4,t_VEC);                  f = carreparfait(ZX_disc(p3)); avma = av;
                 y[2]=lneg(gun); y[1]=lstoi(f? 6: 12);                  return f? _res(6,-1,1): _res(12,-1,3);
               }                }
               else                else
               {                {
                 f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);                  f = carreparfait(ZX_disc(x)); avma = av;
                 if (f) { y[2]=un; y[1]=lstoi(12); }                  return f? _res(12,1,4): _res(24,-1,8);
                 else { y[2]=lneg(gun); y[1]=lstoi(24); }  
               }                }
               y[3]=un; return y;              case 4: avma = av; return _res(6,-1,2);
             case 4: avma=av; y=cgetg(4,t_VEC);  
               y[1]=lstoi(6); y[2]=lneg(gun); y[3]=deux; return y;  
             default: err(bugparier,"galois (bug3)");              default: err(bugparier,"galois (bug3)");
           }            }
         }          }
   
       case 7: z = cgetg(36,t_VEC);        case 7: z = cgetg(36,t_VEC);
           prec = DEFAULTPREC + (long)((gexpo(cb)*7. / BITS_IN_LONG));
         for(;;)          for(;;)
         {          {
           ind = 0; p1=roots(x,prec); p4=gun;            ind = 0; p1=roots(x,prec);
           for (i=1; i<=5; i++)            for (i=1; i<=5; i++)
             for (j=i+1; j<=6; j++)              for (j=i+1; j<=6; j++)
             {              {
               p6 = gadd((GEN)p1[i],(GEN)p1[j]);                GEN t = gadd((GEN)p1[i],(GEN)p1[j]);
               for (k=j+1; k<=7; k++)                for (k=j+1; k<=7; k++) z[++ind] = ladd(t, (GEN)p1[k]);
                 z[++ind] = (long) gadd(p6,(GEN)p1[k]);  
             }              }
           p4 = roots_to_pol(z,0);            p5 = roots_to_ZX(z, &e);
           p5 = grndtoi(greal(p4),&e);  
           e1 = gexpo(gimag(p4)); if (e1>e) e=e1;  
           if (e <= -10) break;            if (e <= -10) break;
           prec = (prec<<1)-2;            prec = (prec<<1)-2;
         }          }
         p6 = modulargcd(p5,derivpol(p5));          if (!ZX_is_squarefree(p5)) goto tchi;
         if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;  
         p2=(GEN)factor(p5)[1];          p2=(GEN)factor(p5)[1];
         switch(lg(p2)-1)          switch(lg(p2)-1)
         {          {
           case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC); y[3]=un;            case 1: f = carreparfait(ZX_disc(x)); avma = av;
             if (f) { y[2]=un; y[1]=lstoi(2520); }              return f? _res(2520,1,6): _res(5040,-1,7);
             else { y[2]=lneg(gun); y[1]=lstoi(5040); }            case 2: f = degpol(p2[1]); avma = av;
             return y;              return (f==7 || f==28)? _res(168,1,5): _res(42,-1,4);
           case 2: f=degpol(p2[1]); avma=av; y=cgetg(4,t_VEC); y[3]=un;            case 3: avma = av; return _res(21,1,3);
             if (f==7 || f==28) { y[2]=un; y[1]=lstoi(168); }            case 4: avma = av; return _res(14,-1,2);
             else { y[2]=lneg(gun); y[1]=lstoi(42); }            case 5: avma = av; return _res(7,1,1);
             return y;  
           case 3: avma=av; y=cgetg(4,t_VEC);  
             y[3]=y[2]=un; y[1]=lstoi(21); return y;  
           case 4: avma=av; y=cgetg(4,t_VEC);  
             y[3]=un; y[2]=lneg(gun); y[1]=lstoi(14); return y;  
           case 5: avma=av; y=cgetg(4,t_VEC);  
             y[3]=y[2]=un; y[1]=lstoi(7); return y;  
           default: err(talker,"galois (bug2)");            default: err(talker,"galois (bug2)");
         }          }
     }      }
     tchi: avma=av1; x=tschirnhaus(x1);      tchi: avma = av1; x = tschirnhaus(x1);
   }    }
 }  }
   
 GEN  GEN
 galoisapply(GEN nf, GEN aut, GEN x)  galoisapply(GEN nf, GEN aut, GEN x)
 {  {
   long av=avma,tetpil,lx,j,N;    gpmem_t av=avma,tetpil;
     long lx,j,N;
   GEN p,p1,y,pol;    GEN p,p1,y,pol;
   
   nf=checknf(nf); pol=(GEN)nf[1];    nf=checknf(nf); pol=(GEN)nf[1];
Line 713  GEN pol_to_monic(GEN pol, GEN *lead);
Line 698  GEN pol_to_monic(GEN pol, GEN *lead);
 int cmp_pol(GEN x, GEN y);  int cmp_pol(GEN x, GEN y);
   
 GEN  GEN
   get_bnfpol(GEN x, GEN *bnf, GEN *nf)
   {
     *bnf = _checkbnf(x);
     *nf  = _checknf(x);
     if (*nf) return (GEN)(*nf)[1];
     if (typ(x) != t_POL) err(idealer1);
     return x;
   }
   
   GEN
 get_nfpol(GEN x, GEN *nf)  get_nfpol(GEN x, GEN *nf)
 {  {
   if (typ(x) == t_POL) { *nf = NULL; return x; }    if (typ(x) == t_POL) { *nf = NULL; return x; }
Line 723  get_nfpol(GEN x, GEN *nf)
Line 718  get_nfpol(GEN x, GEN *nf)
 static GEN  static GEN
 nfiso0(GEN a, GEN b, long fliso)  nfiso0(GEN a, GEN b, long fliso)
 {  {
   ulong av = avma;    gpmem_t av = avma;
   long n,m,i,vb,lx;    long n,m,i,vb,lx;
   GEN nfa,nfb,p1,y,la,lb;    GEN nfa,nfb,p1,y,la,lb;
   
Line 751  nfiso0(GEN a, GEN b, long fliso)
Line 746  nfiso0(GEN a, GEN b, long fliso)
   }    }
   else    else
   {    {
     GEN da = nfa? (GEN)nfa[3]: discsr(a);      GEN da = nfa? (GEN)nfa[3]: ZX_disc(a);
     GEN db = nfb? (GEN)nfb[3]: discsr(b);      GEN db = nfb? (GEN)nfb[3]: ZX_disc(b);
     if (fliso)      if (fliso)
     {      {
       p1=gdiv(da,db);        p1=gdiv(da,db);
Line 787  nfiso0(GEN a, GEN b, long fliso)
Line 782  nfiso0(GEN a, GEN b, long fliso)
     }      }
     y = gen_sort(y, 0, cmp_pol);      y = gen_sort(y, 0, cmp_pol);
     settyp(y, t_VEC);      settyp(y, t_VEC);
     if (vb == 0) delete_var();      if (vb == 0) (void)delete_var();
   }    }
   lx = lg(y); if (lx==1) { avma=av; return gzero; }    lx = lg(y); if (lx==1) { avma=av; return gzero; }
   for (i=1; i<lx; i++)    for (i=1; i<lx; i++)
Line 801  nfiso0(GEN a, GEN b, long fliso)
Line 796  nfiso0(GEN a, GEN b, long fliso)
 }  }
   
 GEN  GEN
 nfisisom(GEN a, GEN b)  nfisisom(GEN a, GEN b) { return nfiso0(a,b,1); }
 {  
   return nfiso0(a,b,1);  
 }  
   
 GEN  GEN
 nfisincl(GEN a, GEN b)  nfisincl(GEN a, GEN b) { return nfiso0(a,b,0); }
 {  
   return nfiso0(a,b,0);  
 }  
   
 /*************************************************************************/  /*************************************************************************/
 /**                                                                     **/  /**                                                                     **/
 /**                            INITALG                                  **/  /**                            INITALG                                  **/
 /**                                                                     **/  /**                                                                     **/
 /*************************************************************************/  /*************************************************************************/
 extern GEN LLL_nfbasis(GEN *x, GEN polr, GEN base, long prec);  
 extern GEN eleval(GEN f,GEN h,GEN a);  
 extern int canon_pol(GEN z);  
 extern GEN mat_to_vecpol(GEN x, long v);  
 extern GEN vecpol_to_mat(GEN v, long n);  
   
   GEN
   get_roots(GEN x,long r1,long prec)
   {
     GEN roo = (typ(x)==t_VEC)? dummycopy(x): roots(x,prec);
     long i, ru = (lg(roo)-1 + r1) >> 1;
   
     for (i=1; i<=r1; i++) roo[i]=lreal((GEN)roo[i]);
     for (   ; i<=ru; i++) roo[i]=roo[(i<<1)-r1];
     roo[0]=evaltyp(t_VEC)|evallg(ru+1); return roo;
   }
   
 /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */  /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
 GEN  GEN
 quicktrace(GEN x, GEN sym)  quicktrace(GEN x, GEN sym)
Line 843  quicktrace(GEN x, GEN sym)
Line 838  quicktrace(GEN x, GEN sym)
 static GEN  static GEN
 trace_col(GEN x, GEN T)  trace_col(GEN x, GEN T)
 {  {
   ulong av = avma;    gpmem_t av = avma;
   GEN t = gzero;    GEN t = gzero;
   long i, l = lg(x);    long i, l = lg(x);
   
Line 853  trace_col(GEN x, GEN T)
Line 848  trace_col(GEN x, GEN T)
   return gerepileuptoint(av, t);    return gerepileuptoint(av, t);
 }  }
   
 /* Seek a new, simpler, polynomial pol defining the same number field as  
  * x (assumed to be monic at this point).  
  * *ptx   receives pol  
  * *ptdx  receives disc(pol)  
  * *ptrev expresses the new root in terms of the old one.  
  * base updated in place.  
  * prec = 0 <==> field is totally real  
  */  
 static void  
 nfinit_reduce(long flag, GEN *px, GEN *pdx, GEN *prev, GEN *pbase, long prec)  
 {  
   GEN chbas,a,phimax,dxn,s,sn,p1,p2,p3,polmax,rev,polr;  
   GEN x = *px, dx = *pdx, base = *pbase;  
   long i,j,nmax,numb,flc,v=varn(x), n=lg(base)-1;  
   
   if (n == 1)  
   {  
     *px = gsub(polx[v],gun); *pdx = gun;  
     *prev = polymodrecip(gmodulcp(gun, x)); return;  
   }  
   polr = prec? roots(x, prec): NULL;  
   if (polr)  
     for (s=gzero,i=1; i<=n; i++) s = gadd(s,gnorm((GEN)polr[i]));  
   else  
     s = subii(sqri((GEN)x[n+1]), shifti((GEN)x[n],1));  
   
   chbas = LLL_nfbasis(&x,polr,base,prec);  
   if (DEBUGLEVEL) msgtimer("LLL basis");  
   
   phimax=polx[v]; polmax=dummycopy(x);  
   nmax=(flag & nf_PARTIAL)?min(n,3):n;  
   
   for (numb=0,i=2; i<=nmax || (!numb && i<=n); i++)  
   { /* cf pols_for_polred */  
     if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); }  
     p1 = a = gmul(base,(GEN)chbas[i]); p3=content(p1);  
     if (gcmp1(p3)) p3 = NULL; else p1 = gdiv(p1,p3);  
     p1 = caract2(x,p1,v);  
     if (p3)  
       for (p2=gun, j=lgef(p1)-2; j>=2; j--)  
       {  
         p2 = gmul(p2,p3); p1[j] = lmul((GEN)p1[j], p2);  
       }  
     p2 = modulargcd(derivpol(p1),p1);  
     if (lgef(p2) > 3) continue;  
   
     if (DEBUGLEVEL>3) outerr(p1);  
     dxn=discsr(p1); flc=absi_cmp(dxn,dx); numb++;  
     if (flc>0) continue;  
   
     if (polr)  
       for (sn=gzero,j=1; j<=n; j++)  
         sn = gadd(sn,gnorm(poleval(a,(GEN)polr[j])));  
     else  
       sn = subii(sqri((GEN)p1[n+1]), shifti((GEN)p1[n],1));  
     if (flc>=0)  
     {  
       flc=gcmp(sn,s);  
       if (flc>0 || (!flc && gpolcomp(p1,polmax) >= 0)) continue;  
     }  
     dx=dxn; s=sn; polmax=p1; phimax=a;  
   }  
   if (!numb)  
   {  
     if (gisirreducible(x)!=gun) err(redpoler,"nfinit_reduce");  
     err(talker,"you have found a counter-example to a conjecture, "  
                "please send us\nthe polynomial as soon as possible");  
   }  
   if (phimax == polx[v]) /* no improvement */  
     rev = gmodulcp(phimax,x);  
   else  
   {  
     if (canon_pol(polmax) < 0) phimax = gneg_i(phimax);  
     if (DEBUGLEVEL>1) fprintferr("polmax = %Z\n",polmax);  
     p2 = gmodulcp(phimax,x); rev = polymodrecip(p2);  
     a = base; base = cgetg(n+1,t_VEC);  
     p1 = (GEN)rev[2];  
     for (i=1; i<=n; i++)  
       base[i] = (long)eleval(polmax, (GEN)a[i], p1);  
     p1 = vecpol_to_mat(base,n); p2=denom(p1); p1=gmul(p2,p1);  
     p1 = gdiv(hnfmodid(p1,p2), p2);  
     base = mat_to_vecpol(p1,v);  
   }  
   *px=polmax; *pdx=dx; *prev=rev; *pbase = base;  
 }  
   
 /* pol belonging to Z[x], return a monic polynomial generating the same field  /* pol belonging to Z[x], return a monic polynomial generating the same field
  * as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing   * as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing
  * by the content), and to to leading coeff otherwise.   * by the content), and to to leading coeff otherwise.
Line 948  GEN
Line 857  GEN
 pol_to_monic(GEN pol, GEN *lead)  pol_to_monic(GEN pol, GEN *lead)
 {  {
   long n = lgef(pol)-1;    long n = lgef(pol)-1;
   GEN p1;  
   
   if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; }    if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; }
     return primitive_pol_to_monic(primpart(pol), lead);
   p1=content(pol); if (!gcmp1(p1)) pol = gdiv(pol,p1);  
   return primitive_pol_to_monic(pol,lead);  
 }  }
   
 /* NB: TI integral, det(TI) = d_K, ideal/dideal = codifferent */  /* Let T = (Tr(wi wj)), then T^(-1) Z^n = codifferent = coD/dcoD
 GEN   * NB: det(T) = dK, TI = dK T^(-1) integral */
 make_MDI(GEN nf, GEN TI, GEN *ideal, GEN *dideal)  static GEN
   make_MDI(GEN nf, GEN TI, GEN *coD, GEN *dcoD)
 {  {
   GEN c = content(TI);    GEN c, MDI, dK = (GEN)nf[3];
   GEN p1;  
   
   *dideal = divii((GEN)nf[3], c);    TI = Q_primitive_part(TI, &c);
   *ideal = p1 = hnfmodid(gdiv(TI,c), *dideal);    *dcoD = c? diviiexact(dK, c): dK;
   return gmul(c, ideal_two_elt(nf, p1));    *coD = hnfmodid(TI, *dcoD);
     MDI = ideal_two_elt(nf, *coD);
     return c? gmul(c, MDI): MDI;
 }  }
   
 /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */  
 GEN  
 make_M(GEN basden,GEN roo)  
 {  
   GEN p1,d,M, bas=(GEN)basden[1], den=(GEN)basden[2];  
   long i,j, ru = lg(roo), n = lg(bas);  
   M = cgetg(n,t_MAT);  
   for (j=1; j<n; j++)  
   {  
     p1=cgetg(ru,t_COL); M[j]=(long)p1;  
     for (i=1; i<ru; i++)  
       p1[i]=(long)poleval((GEN)bas[j],(GEN)roo[i]);  
   }  
   if (den)  
   {  
     long prec = precision((GEN)roo[1]);  
     GEN invd,rd = cgetr(prec+1);  
     for (j=1; j<n; j++)  
     {  
       d = (GEN)den[j]; if (!d) continue;  
       p1 = (GEN)M[j]; affir(d,rd); invd = ginv(rd);  
       for (i=1; i<ru; i++) p1[i] = lmul((GEN)p1[i], invd);  
     }  
   }  
   if (DEBUGLEVEL>4) msgtimer("matrix M");  
   return M;  
 }  
   
 GEN  
 make_MC(long r1,GEN M)  
 {  
   long i,j,av,tetpil, n = lg(M), ru = lg(M[1]);  
   GEN p1,p2,MC=cgetg(ru,t_MAT);  
   
   for (j=1; j<ru; j++)  
   {  
     p1=cgetg(n,t_COL); MC[j]=(long)p1;  
     for (i=1; i<n; i++)  
     {  
       av=avma; p2=gconj(gcoeff(M,j,i)); tetpil=avma;  
       p1[i] = (j<=r1)? (long)p2: lpile(av,tetpil,gmul2n(p2,1));  
     }  
   }  
   if (DEBUGLEVEL>4) msgtimer("matrix MC");  
   return MC;  
 }  
   
 GEN  
 get_roots(GEN x,long r1,long ru,long prec)  
 {  
   GEN roo = (typ(x)==t_VEC)? dummycopy(x): roots(x,prec);  
   long i;  
   
   for (i=1; i<=r1; i++) roo[i]=lreal((GEN)roo[i]);  
   for (   ; i<=ru; i++) roo[i]=roo[(i<<1)-r1];  
   roo[0]=evaltyp(t_VEC)|evallg(ru+1); return roo;  
 }  
   
 /* assume lg(nf) > 3 && typ(nf) = container [hopefully a genuine nf] */  /* assume lg(nf) > 3 && typ(nf) = container [hopefully a genuine nf] */
 long  long
 nf_get_r1(GEN nf)  nf_get_r1(GEN nf)
Line 1043  nf_get_r2(GEN nf)
Line 893  nf_get_r2(GEN nf)
     err(talker,"false nf in nf_get_r2");      err(talker,"false nf in nf_get_r2");
   return itos((GEN)x[2]);    return itos((GEN)x[2]);
 }  }
   void
   nf_get_sign(GEN nf, long *r1, long *r2)
   {
     GEN x = (GEN)nf[2];
     if (typ(x) != t_VEC || lg(x) != 3
         || typ(x[1]) != t_INT || typ(x[2]) != t_INT)
       err(talker,"false nf in nf_get_sign");
     *r1 = itos((GEN)x[1]);
     *r2 = (degpol(nf[1]) - *r1) >> 1;
   }
   
 extern GEN mulmat_pol(GEN A, GEN x);  
   
 static GEN  static GEN
 get_T(GEN mul, GEN x, GEN bas, GEN den)  get_Tr(GEN mul, GEN x, GEN basden)
 {  {
     GEN tr,T,T1,sym, bas = (GEN)basden[1], den = (GEN)basden[2];
   long i,j,n = lg(bas)-1;    long i,j,n = lg(bas)-1;
   GEN tr,T,T1,sym;  
   T = cgetg(n+1,t_MAT);    T = cgetg(n+1,t_MAT);
   T1 = cgetg(n+1,t_COL);    T1 = cgetg(n+1,t_COL);
   sym = polsym(x, n-1);    sym = polsym(x, n-1);
Line 1076  get_T(GEN mul, GEN x, GEN bas, GEN den)
Line 934  get_T(GEN mul, GEN x, GEN bas, GEN den)
 GEN  GEN
 get_bas_den(GEN bas)  get_bas_den(GEN bas)
 {  {
   GEN z,d,den, dbas = dummycopy(bas);    GEN z,b,d,den, dbas = dummycopy(bas);
   long i, c = 0, l = lg(bas);    long i, l = lg(bas);
     int power = 1;
   den = cgetg(l,t_VEC);    den = cgetg(l,t_VEC);
   for (i=1; i<l; i++)    for (i=1; i<l; i++)
   {    {
     d = denom(content((GEN)dbas[i]));      b = Q_remove_denom((GEN)bas[i], &d);
     if (is_pm1(d)) d = NULL; else { dbas[i] = lmul((GEN)dbas[i],d); c++; }      dbas[i]= (long)b;
     den[i] = (long)d;      den[i] = (long)d; if (d) power = 0;
   }    }
   if (!c) den = NULL; /* power basis */    if (power) den = NULL; /* power basis */
   z = cgetg(3,t_VEC);    z = cgetg(3,t_VEC);
   z[1] = (long)dbas;    z[1] = (long)dbas;
   z[2] = (long)den; return z;    z[2] = (long)den; return z;
Line 1101  _mulii(GEN x, GEN y)
Line 960  _mulii(GEN x, GEN y)
 }  }
   
 GEN  GEN
 get_mul_table(GEN x,GEN basden,GEN invbas,GEN *T)  get_mul_table(GEN x,GEN basden,GEN invbas)
 {  {
   long i,j, n = degpol(x);    long i,j, n = degpol(x);
   GEN z,d, mul = cgetg(n*n+1,t_MAT), bas=(GEN)basden[1], den=(GEN)basden[2];    GEN z,d, mul = cgetg(n*n+1,t_MAT), bas=(GEN)basden[1], den=(GEN)basden[2];
   
   for (j=1; j<=n*n; j++) mul[j]=lgetg(n+1,t_COL);    for (j=1; j<=n*n; j++) mul[j]=lgetg(n+1,t_COL);
   for (i=1; i<=n; i++)    for (i=1; i<=n; i++)
     for (j=i; j<=n; j++)      for (j=i; j<=n; j++)
     {      {
       ulong av = avma;        gpmem_t av = avma;
       z = gres(gmul((GEN)bas[j],(GEN)bas[i]), x);        z = gres(gmul((GEN)bas[j],(GEN)bas[i]), x);
       z = mulmat_pol(invbas, z); /* integral column */        z = mulmat_pol(invbas, z); /* integral column */
       if (den)        if (den)
Line 1120  get_mul_table(GEN x,GEN basden,GEN invbas,GEN *T)
Line 979  get_mul_table(GEN x,GEN basden,GEN invbas,GEN *T)
       }        }
       mul[j+(i-1)*n] = mul[i+(j-1)*n] = lpileupto(av,z);        mul[j+(i-1)*n] = mul[i+(j-1)*n] = lpileupto(av,z);
     }      }
   if (T) *T = get_T(mul,x,bas,den);  
   return mul;    return mul;
 }  }
   
 /* fill mat = nf[5], as well as nf[8] and nf[9]  /* as get_Tr, mul_table not precomputed */
  * If (small) only compute a subset (use dummy 0s for the rest) */  static GEN
   make_Tr(GEN x, GEN w)
   {
     long i,j, n = degpol(x);
     GEN p1,p2,t,d;
     GEN sym = cgetg(n+2,t_VEC);
     GEN den = cgetg(n+1,t_VEC);
     GEN T = cgetg(n+1,t_MAT);
     gpmem_t av;
   
     sym = polsym(x, n-1);
     p1 = get_bas_den(w);
     w   = (GEN)p1[1];
     den = (GEN)p1[2];
     for (i=1; i<=n; i++)
     {
       p1 = cgetg(n+1,t_COL); T[i] = (long)p1;
       for (j=1; j<i ; j++) p1[j] = coeff(T,i,j);
       for (   ; j<=n; j++)
       {
         av = avma;
         p2 = gres(gmul((GEN)w[i],(GEN)w[j]), x);
         t = quicktrace(p2, sym);
         if (den)
         {
           d = _mulii((GEN)den[i],(GEN)den[j]);
           if (d) t = diviiexact(t, d);
         }
         p1[j] = lpileuptoint(av, t);
       }
     }
     return T;
   }
   
   /* compute roots so that _absolute_ accuracy of M >= prec [also holds for G] */
   static void
   get_roots_for_M(nffp_t *F)
   {
     long n, eBD, er, PREC;
   
     if (F->extraprec < 0)
     { /* not initialized yet */
       n = degpol(F->x);
       eBD = 1 + gexpo((GEN)F->basden[1]);
       er  = 1 + gexpo(F->ro? F->ro: cauchy_bound(F->x));
       if (er < 0) er = 0;
       F->extraprec = ((n*er + eBD + (long)log2(n)) >> TWOPOTBITS_IN_LONG);
     }
   
     PREC = F->prec + F->extraprec;
     if (F->ro && gprecision((GEN)F->ro[1]) >= PREC) return;
     F->ro = get_roots(F->x, F->r1, PREC);
   }
   
   /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
   static void
   make_M(nffp_t *F, int trunc)
   {
     GEN bas = (GEN)F->basden[1], den = (GEN)F->basden[2], ro = F->ro;
     GEN m, d, M;
     long i, j, l = lg(ro), n = lg(bas);
     M = cgetg(n,t_MAT);
       m = cgetg(l, t_COL); M[1] = (long)m;
       for (i=1; i<l; i++) m[i] = un; /* bas[1] = 1 */
     for (j=2; j<n; j++)
     {
       m = cgetg(l,t_COL); M[j] = (long)m;
       for (i=1; i<l; i++) m[i] = (long)poleval((GEN)bas[j], (GEN)ro[i]);
     }
     if (den)
     {
       GEN invd, rd = cgetr(F->prec + F->extraprec);
       for (j=2; j<n; j++)
       {
         d = (GEN)den[j]; if (!d) continue;
         m = (GEN)M[j]; affir(d,rd); invd = ginv(rd);
         for (i=1; i<l; i++) m[i] = lmul((GEN)m[i], invd);
       }
     }
   
     if (trunc && gprecision(M) > F->prec)
     {
       M     = gprec_w(M, F->prec);
       F->ro = gprec_w(ro,F->prec);
     }
     if (DEBUGLEVEL>4) msgtimer("matrix M");
     F->M = M;
   }
   
   /* return G real such that G~ * G = T_2 */
   static void
   make_G(nffp_t *F)
   {
     GEN G, m, g, r, M = F->M, sqrt2 = gsqrt(gdeux, F->prec + F->extraprec);
     long i, j, k, r1 = F->r1, l = lg(M);
   
     G = cgetg(l, t_MAT);
     for (j=1; j<l; j++)
     {
       g = cgetg(l, t_COL);
       G[j] = (long)g; m = (GEN)M[j];
       for (k=i=1; i<=r1; i++) g[k++] = m[i];
       for (     ; k < l; i++)
       {
         r = (GEN)m[i];
         if (typ(r) == t_COMPLEX)
         {
           g[k++] = lmpmul(sqrt2, (GEN)r[1]);
           g[k++] = lmpmul(sqrt2, (GEN)r[2]);
         }
         else
         {
           g[k++] = lmpmul(sqrt2, r);
           g[k++] = zero;
         }
       }
     }
     F->G = G;
   }
   
   static void
   make_M_G(nffp_t *F, int trunc)
   {
     get_roots_for_M(F);
     make_M(F, trunc);
     make_G(F);
   }
   
 void  void
 get_nf_matrices(GEN nf, long small)  remake_GM(GEN nf, nffp_t *F, long prec)
 {  {
   GEN x=(GEN)nf[1],dK=(GEN)nf[3],index=(GEN)nf[4],ro=(GEN)nf[6],bas=(GEN)nf[7];    F->x  = (GEN)nf[1];
   GEN basden,mul,invbas,M,MC,T,MDI,D,TI,A,dA,mat;    F->ro = NULL;
   long r1 = nf_get_r1(nf), n = lg(bas)-1;    F->r1 = nf_get_r1(nf);
     F->basden = get_bas_den((GEN)nf[7]);
     F->extraprec = -1;
     F->prec = prec; make_M_G(F, 1);
   }
   
   mat = cgetg(small? 4: 8,t_VEC); nf[5] = (long)mat;  static void
   basden = get_bas_den(bas);  nffp_init(nffp_t *F, nfbasic_t *T, GEN ro, long prec)
   M = make_M(basden,ro);  {
   MC = make_MC(r1,M);    F->x  = T->x;
   mat[1]=(long)M;    F->ro = ro;
   mat[3]=(long)mulmat_real(MC,M);    F->r1 = T->r1;
   if (small) { nf[8]=nf[9]=mat[2]=zero; return; }    if (!T->basden) T->basden = get_bas_den(T->bas);
     F->basden = T->basden;
     F->extraprec = -1;
     F->prec = prec;
   }
   
   invbas = QM_inv(vecpol_to_mat(bas,n), gun);  static void
   mul = get_mul_table(x,basden,invbas,&T);  get_nf_fp_compo(nfbasic_t *T, nffp_t *F, GEN ro, long prec)
   {
     nffp_init(F,T,ro,prec);
     make_M_G(F, 1);
   }
   
   static GEN
   get_sign(long r1, long n)
   {
     GEN s = cgetg(3, t_VEC);
     s[1] = lstoi(r1);
     s[2] = lstoi((n - r1) >> 1); return s;
   }
   
   GEN
   nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec)
   {
     GEN nf = cgetg(10,t_VEC);
     GEN x = T->x;
     GEN invbas,Tr,D,TI,A,dA, mat = cgetg(8,t_VEC);
     nffp_t F;
     get_nf_fp_compo(T, &F, ro, prec);
   
     nf[1] = (long)T->x;
     nf[2] = (long)get_sign(T->r1, degpol(T->x));
     nf[3] = (long)T->dK;
     nf[4] = (long)T->index;
     nf[6] = (long)F.ro;
     nf[5] = (long)mat;
     nf[7] = (long)T->bas;
   
     mat[1] = (long)F.M;
     mat[2] = (long)F.G;
   
     invbas = QM_inv(vecpol_to_mat(T->bas, lg(T->bas)-1), gun);
     nf[8] = (long)invbas;
     nf[9] = (long)get_mul_table(x, F.basden, invbas);
   if (DEBUGLEVEL) msgtimer("mult. table");    if (DEBUGLEVEL) msgtimer("mult. table");
   
   nf[8]=(long)invbas;  
   nf[9]=(long)mul;  
   
   TI = ZM_inv(T, dK); /* dK T^-1 */    Tr = get_Tr((GEN)nf[9], x, F.basden);
   MDI = make_MDI(nf,TI, &A, &dA);    TI = ZM_inv(Tr, T->dK); /* dK T^-1 */
   mat[6]=(long)TI;    mat[6] = (long)TI;
   mat[7]=(long)MDI; /* needed in idealinv below */    mat[7] = (long)make_MDI(nf,TI, &A, &dA); /* needed in idealinv below */
   if (is_pm1(index))    if (is_pm1(T->index))
     D = idealhermite_aux(nf, derivpol(x));      D = idealhermite_aux(nf, derivpol(x));
   else    else
     D = gmul(dA, idealinv(nf, A));      D = gmul(dA, idealinv(nf, A));
   mat[2]=(long)MC;    mat[3] = zero; /* FIXME: was gram matrix of current mat[2]. Useless */
   mat[4]=(long)T;    mat[4] = (long)Tr;
   mat[5]=(long)D;    mat[5] = (long)D;
   if (DEBUGLEVEL) msgtimer("matrices");    if (DEBUGLEVEL) msgtimer("matrices");
     return nf;
 }  }
   
 /* Initialize the number field defined by the polynomial x (in variable v)  static GEN
  * flag & nf_REGULAR  hnffromLLL(GEN nf)
  *    regular behaviour.  {
  * flag & nf_SMALL    GEN d, x;
  *    compute only nf[1] (pol), nf[2] (signature), nf[5][3] (T2) and    x = vecpol_to_mat((GEN)nf[7], degpol(nf[1]));
  *    nf[7] (integer basis), the other components are filled with gzero.    x = Q_remove_denom(x, &d);
  * flag & nf_REDUCE    if (!d) return x; /* power basis */
  *    try a polred first.    return gauss(hnfmodid(x, d), x);
  * flag & nf_PARTIAL  }
  *    do a partial polred, not a polredabs  
  * flag & nf_ORIG  
  *    do a polred and return [nfinit(x),Mod(a,red)], where  
  *    Mod(a,red)=Mod(v,x) (i.e return the base change).  
  */  
   
 /* here x can be a polynomial, an nf or a bnf */  static GEN
   nfbasechange(GEN u, GEN x)
   {
     long i,lx;
     GEN y;
     switch(typ(x))
     {
       case t_COL: /* nfelt */
         return gmul(u, x);
   
       case t_MAT: /* ideal */
         y = dummycopy(x);
         lx = lg(x);
         for (i=1; i<lx; i++) y[i] = lmul(u, (GEN)y[i]);
         break;
   
       case t_VEC: /* pr */
         checkprimeid(x);
         y = dummycopy(x);
         y[2] = lmul(u, (GEN)y[2]);
         y[5] = lmul(u, (GEN)y[5]);
         break;
       default: y = x;
     }
     return y;
   }
   
 GEN  GEN
 initalgall0(GEN x, long flag, long prec)  nffromhnfbasis(GEN nf, GEN x)
 {  {
   GEN lead = NULL,nf,ro,bas,mat,rev,dK,dx,index,res;    long tx = typ(x);
   long av=avma,n,i,r1,r2,ru,PRECREG;    gpmem_t av = avma;
     GEN u;
     if (!is_vec_t(tx)) return gcopy(x);
     nf = checknf(nf);
     u = hnffromLLL(nf);
     return gerepilecopy(av, nfbasechange(u, x));
   }
   
   if (DEBUGLEVEL) timer2();  GEN
   if (typ(x)==t_POL)  nftohnfbasis(GEN nf, GEN x)
   {
     long tx = typ(x);
     gpmem_t av = avma;
     GEN u;
     if (!is_vec_t(tx)) return gcopy(x);
     nf = checknf(nf);
     u = ZM_inv(hnffromLLL(nf), gun);
     return gerepilecopy(av, nfbasechange(u, x));
   }
   
   static GEN
   get_red_G(nfbasic_t *T, GEN *pro)
   {
     GEN G, u, u0 = NULL;
     gpmem_t av;
     long i, prec, extraprec, n = degpol(T->x);
     nffp_t F;
   
     extraprec = (long) (0.5 * n * sizeof(long) / 8.);
     prec = DEFAULTPREC + extraprec;
     nffp_init(&F, T, *pro, prec);
     av = avma;
     for (i=1; ; i++)
   {    {
     n=degpol(x); if (n<=0) err(constpoler,"initalgall0");      F.prec = prec; make_M_G(&F, 0); G = F.G;
     check_pol_int(x,"initalg");      if (u0) G = gmul(G, u0);
     if (gisirreducible(x) == gzero) err(redpoler,"nfinit");      if (DEBUGLEVEL)
     if (!gcmp1((GEN)x[n+2]))        fprintferr("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
                     prec + F.extraprec, prec, F.extraprec);
       if ((u = lllfp_marked(1, G, 100, 2, prec, 0)))
     {      {
       x = pol_to_monic(x,&lead);        if (typ(u) == t_MAT) break;
       if (!(flag & nf_SMALL))        u = (GEN)u[1];
       {        if (u0) u0 = gerepileupto(av, gmul(u0,u));
         if (!(flag & nf_REDUCE))        else    u0 = gerepilecopy(av, u);
           err(warner,"non-monic polynomial. Result of the form [nf,c]");  
         flag = flag | nf_REDUCE | nf_ORIG;  
       }  
     }      }
     bas = allbase4(x,0,&dK,NULL);      if (i == MAXITERPOL) err(accurer,"red_T2");
     if (DEBUGLEVEL) msgtimer("round4");      prec = (prec<<1)-2 + (gexpo(u0) >> TWOPOTBITS_IN_LONG);
     dx = discsr(x); r1 = sturm(x);      F.ro = NULL;
       if (DEBUGLEVEL) err(warnprec,"get_red_G", prec);
   }    }
     *pro = F.ro; return u0? gmul(u0,u): u;
   }
   
   /* Return the base change matrix giving an LLL-reduced basis for the
    * integer basis of nf(T).
    * *ro = roots of x, computed to precision prec [or NULL] */
   static GEN
   get_LLL_basis(nfbasic_t *T, GEN *pro)
   {
     GEN u;
     if (T->r1 != degpol(T->x)) u = get_red_G(T, pro);
   else    else
   {    {
     i = lg(x);      u = lllint_marked(1, make_Tr(T->x, T->bas), 100, 1, &u,NULL,NULL);
     if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)      if (!u) u = idmat(1);
     { /* polynomial + integer basis */  
       bas=(GEN)x[2]; x=(GEN)x[1]; n=degpol(x);  
       if (typ(bas) == t_MAT) { mat = bas; bas = mat_to_vecpol(mat,varn(x)); }  
       else mat = vecpol_to_mat(bas,n);  
       dx = discsr(x); r1 = sturm(x);  
       dK = gmul(dx, gsqr(det2(mat)));  
     }  
     else  
     {  
       nf=checknf(x);  
       bas=(GEN)nf[7]; x=(GEN)nf[1]; n=degpol(x);  
       dK=(GEN)nf[3]; dx=mulii(dK, sqri((GEN)nf[4]));  
       r1 = nf_get_r1(nf);  
     }  
     bas[1]=lpolun[varn(x)]; /* it may be gun => SEGV later */  
   }    }
   r2=(n-r1)>>1; ru=r1+r2;    return u;
   PRECREG = prec + (expi(dK)>>(TWOPOTBITS_IN_LONG+1))  }
                  + (long)((sqrt((double)n)+3) / sizeof(long) * 4);  
   rev = NULL;  static void
   if (flag & nf_REDUCE)  set_LLL_basis(nfbasic_t *T, GEN *pro)
   {
     T->bas = gmul(T->bas, get_LLL_basis(T, pro));
     T->basden = NULL; /* recompute */
   }
   
   typedef struct {
     GEN xbest, dxbest;
     long ind, indmax, indbest;
   } ok_pol_t;
   
   /* is xn better than x ? */
   static int
   better_pol(GEN xn, GEN dxn, GEN x, GEN dx)
   {
     int fl;
     if (!x) return 1;
     fl = absi_cmp(dxn, dx);
     return (fl < 0 || (!fl && gpolcomp(xn, x) < 0));
   }
   
   static GEN
   ok_pol(void *TT, GEN xn)
   {
     ok_pol_t *T = (ok_pol_t*)TT;
     GEN dxn;
   
     if (++T->ind > T->indmax && T->xbest) return T->xbest;
   
     if (!ZX_is_squarefree(xn)) return (T->ind == T->indmax)? T->xbest: NULL;
     if (DEBUGLEVEL>3) outerr(xn);
     dxn = ZX_disc(xn);
     if (better_pol(xn, dxn, T->xbest, T->dxbest))
   {    {
     nfinit_reduce(flag, &x, &dx, &rev, &bas, r1==n? 0: prec);      T->dxbest = dxn; T->xbest = xn; T->indbest = T->ind;
     if (DEBUGLEVEL) msgtimer("polred");  
   }    }
   if (!carrecomplet(divii(dx,dK),&index))    if (T->ind >= T->indmax) return T->xbest;
     err(bugparier,"nfinit (incorrect discriminant)");    return NULL;
   }
   
   ro=get_roots(x,r1,ru,PRECREG);  /* z in Z[X] with positive leading coeff. Set z := z(-X) or z(X) so that
   if (DEBUGLEVEL) msgtimer("roots");   * second coeff is > 0. In place. */
   static int
   canon_pol(GEN z)
   {
     long i,s;
   
   nf=cgetg(10,t_VEC);    for (i=lgef(z)-2; i>=2; i-=2)
   nf[1]=(long)x;  
   nf[2]=lgetg(3,t_VEC);  
   mael(nf,2,1)=lstoi(r1);  
   mael(nf,2,2)=lstoi(r2);  
   nf[3]=(long)dK;  
   nf[4]=(long)index;  
   nf[6]=(long)ro;  
   nf[7]=(long)bas;  
   get_nf_matrices(nf, flag & nf_SMALL);  
   
   if (flag & nf_ORIG)  
   {    {
     if (!rev) err(talker,"bad flag in initalgall0");      s = signe(z[i]);
     res = cgetg(3,t_VEC);      if (s > 0)
     res[1]=(long)nf; nf = res;      {
     res[2]=lead? ldiv(rev,lead): (long)rev;        for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]);
         return -1;
       }
       if (s) return 1;
   }    }
   return gerepilecopy(av, nf);    return 0;
 }  }
   
   static GEN _polred(GEN x, GEN a, GEN *pta, FP_chk_fun *CHECK);
   
   /* Seek a simpler, polynomial pol defining the same number field as
    * x (assumed to be monic at this point) */
   static GEN
   nfpolred(int part, nfbasic_t *T)
   {
     GEN x = T->x, dx = T->dx, a = T->bas;
     GEN phi, xbest, dxbest, mat, d, rev;
     long i, n = lg(a)-1, v = varn(x);
     ok_pol_t O;
     FP_chk_fun chk;
   
     if (degpol(x) == 1) { T->x = gsub(polx[v],gun); return gun; }
   
     if (!dx) dx = mulii(T->dK, sqri(T->index));
   
     O.ind    = 0;
     O.indmax = part? min(n,3): n;
     O.xbest  = NULL;
     chk.f    = &ok_pol;
     chk.data = (void*)&O;
     if (!_polred(x, a, NULL, &chk))
       err(talker,"you found a counter-example to a conjecture, please report!");
     xbest = O.xbest; dxbest = O.dxbest;
     if (!better_pol(xbest, dxbest, x, dx)) return NULL; /* no improvement */
   
     /* update T */
     phi = (GEN)a[O.indbest];
     if (canon_pol(xbest) < 0) phi = gneg_i(phi);
     if (DEBUGLEVEL>1) fprintferr("xbest = %Z\n",xbest);
     rev = modreverse_i(phi, x);
     for (i=1; i<=n; i++) a[i] = (long)RX_RXQ_compo((GEN)a[i], rev, xbest);
     mat = vecpol_to_mat(Q_remove_denom(a, &d), n);
     if (!is_pm1(d)) mat = gdiv(hnfmodid(mat,d), d); else mat = idmat(n);
   
     (void)carrecomplet(diviiexact(dxbest,T->dK), &(T->index));
     T->bas= mat_to_vecpol(mat, v);
     T->dx = dxbest;
     T->x  = xbest; return rev;
   }
   
 GEN  GEN
 initalgred(GEN x, long prec)  get_nfindex(GEN bas)
 {  {
   return initalgall0(x,nf_REDUCE,prec);    gpmem_t av = avma;
     long n = lg(bas)-1;
     GEN d, mat = vecpol_to_mat(Q_remove_denom(bas, &d), n);
     if (!d) { avma = av; return gun; }
     return gerepileuptoint(av, diviiexact(gpowgs(d, n), det(mat)));
 }  }
   
   void
   nfbasic_init(GEN x, long flag, GEN fa, nfbasic_t *T)
   {
     GEN bas, dK, dx, index;
     long r1;
   
     T->basden = NULL;
     T->lead   = NULL;
     if (DEBUGLEVEL) (void)timer2();
     if (typ(x) == t_POL)
     {
       check_pol_int(x, "nfinit");
       if (gisirreducible(x) == gzero) err(redpoler, "nfinit");
       x = pol_to_monic(x, &(T->lead));
       bas = allbase(x, flag, &dx, &dK, &index, &fa);
       if (DEBUGLEVEL) msgtimer("round4");
       r1 = sturm(x);
     }
     else if (typ(x) == t_VEC && lg(x)==3 && typ(x[1])==t_POL)
     { /* monic integral polynomial + integer basis */
       GEN mat;
       bas = (GEN)x[2]; x = (GEN)x[1];
       if (typ(bas) == t_MAT)
         { mat = bas; bas = mat_to_vecpol(mat,varn(x)); }
       else
           mat = vecpol_to_mat(bas, lg(bas)-1);
       index = get_nfindex(bas);
       dx = ZX_disc(x);
       dK = diviiexact(dx, sqri(index));
       r1 = sturm(x);
     }
     else
     { /* nf, bnf, bnr */
       GEN nf = checknf(x);
       x     = (GEN)nf[1];
       dK    = (GEN)nf[3];
       index = (GEN)nf[4];
       bas   = (GEN)nf[7];
       r1 = nf_get_r1(nf); dx = NULL;
     }
   
     T->x     = x;
     T->r1    = r1;
     T->dx    = dx;
     T->dK    = dK;
     T->bas   = bas;
     T->index = index;
   }
   
   /* Initialize the number field defined by the polynomial x (in variable v)
    * flag & nf_RED:     try a polred first.
    * flag & nf_PARTRED: as nf_RED but check the first two elements only.
    * flag & nf_ORIG
    *    do a polred and return [nfinit(x), Mod(a,red)], where
    *    Mod(a,red) = Mod(v,x) (i.e return the base change). */
 GEN  GEN
 initalgred2(GEN x, long prec)  _initalg(GEN x, long flag, long prec)
 {  {
   return initalgall0(x,nf_REDUCE|nf_ORIG,prec);    const gpmem_t av = avma;
     GEN nf, rev = NULL, ro = NULL;
     nfbasic_t T;
   
     nfbasic_init(x, flag, NULL, &T);
   
     if (T.lead && !(flag & (nf_RED|nf_PARTRED)))
     {
       err(warner,"non-monic polynomial. Result of the form [nf,c]");
       flag |= nf_PARTRED | nf_ORIG;
     }
     if (flag & (nf_RED|nf_PARTRED))
     {
       set_LLL_basis(&T, &ro);
       rev = nfpolred(flag & nf_PARTRED, &T);
       if (rev) ro = NULL; /* changed T.x */
       if (flag & nf_ORIG)
       {
         if (!rev) rev = polx[varn(T.x)]; /* no improvement */
         if (T.lead) rev = gdiv(rev, T.lead);
         rev = to_polmod(rev, T.x);
       }
       if (DEBUGLEVEL) msgtimer("polred");
     }
   
     set_LLL_basis(&T, &ro);
     if (DEBUGLEVEL) msgtimer("LLL basis");
   
     nf = nfbasic_to_nf(&T, ro, prec);
     if (flag & nf_ORIG)
     {
       GEN res = cgetg(3,t_VEC);
       res[1] = (long)nf;
       res[2] = (long)rev; nf = res;
     }
     return gerepilecopy(av, nf);
 }  }
   
 GEN  GEN
   initalgred(GEN x, long prec)  { return _initalg(x, nf_RED, prec); }
   GEN
   initalgred2(GEN x, long prec) { return _initalg(x, nf_RED|nf_ORIG, prec); }
   GEN
   initalg(GEN x, long prec)     { return _initalg(x, 0, prec); }
   
   GEN
 nfinit0(GEN x, long flag,long prec)  nfinit0(GEN x, long flag,long prec)
 {  {
   switch(flag)    switch(flag)
   {    {
     case 0:      case 0:
     case 1: return initalgall0(x,nf_REGULAR,prec);      case 1: return _initalg(x,0,prec);
     case 2: return initalgall0(x,nf_REDUCE,prec);      case 2: return _initalg(x,nf_RED,prec);
     case 3: return initalgall0(x,nf_REDUCE|nf_ORIG,prec);      case 3: return _initalg(x,nf_RED|nf_ORIG,prec);
     case 4: return initalgall0(x,nf_REDUCE|nf_PARTIAL,prec);      case 4: return _initalg(x,nf_PARTRED,prec);
     case 5: return initalgall0(x,nf_REDUCE|nf_ORIG|nf_PARTIAL,prec);      case 5: return _initalg(x,nf_PARTRED|nf_ORIG,prec);
     case 6: return initalgall0(x,nf_SMALL,prec);  
     default: err(flagerr,"nfinit");      default: err(flagerr,"nfinit");
   }    }
   return NULL; /* not reached */    return NULL; /* not reached */
 }  }
   
 GEN  
 initalg(GEN x, long prec)  
 {  
   return initalgall0(x,nf_REGULAR,prec);  
 }  
   
 /* assume x a bnr/bnf/nf */  /* assume x a bnr/bnf/nf */
 long  long
 nfgetprec(GEN x)  nfgetprec(GEN x)
 {  {
   GEN nf = checknf(x), ro = (GEN)nf[6];    GEN nf = checknf(x), ro = (GEN)nf[6];
   return (typ(ro)==t_VEC)?precision((GEN)ro[1]): DEFAULTPREC;    return (typ(ro)==t_VEC)? precision((GEN)ro[1]): DEFAULTPREC;
 }  }
   
 GEN  GEN
 nfnewprec(GEN nf, long prec)  nfnewprec(GEN nf, long prec)
 {  {
   long av=avma,r1,r2,ru,n;    gpmem_t av = avma;
   GEN y,pol,ro,basden,MC,mat,M;    GEN NF;
     nffp_t F;
   
   if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec");    if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec");
   if (lg(nf) == 11) return bnfnewprec(nf,prec);    if (lg(nf) == 11) return bnfnewprec(nf,prec);
   if (lg(nf) ==  7) return bnrnewprec(nf,prec);    if (lg(nf) ==  7) return bnrnewprec(nf,prec);
   (void)checknf(nf);    (void)checknf(nf);
   y = dummycopy(nf);    NF = dummycopy(nf);
   pol=(GEN)nf[1]; n=degpol(pol);    NF[5] = (long)dummycopy((GEN)NF[5]);
   r1 = nf_get_r1(nf);    remake_GM(NF, &F, prec);
   r2 = (n - r1) >> 1; ru = r1+r2;    NF[6]  = (long)F.ro;
   mat=dummycopy((GEN)nf[5]);    mael(NF,5,1) = (long)F.M;
   ro=get_roots(pol,r1,ru,prec);    mael(NF,5,2) = (long)F.G;
   y[5]=(long)mat;    return gerepilecopy(av, NF);
   y[6]=(long)ro;  }
   basden = get_bas_den((GEN)nf[7]);  
   M = make_M(basden,ro);  /********************************************************************/
   MC = make_MC(r1,M);  /**                                                                **/
   mat[1]=(long)M;  /**                           POLRED                               **/
   if (typ(nf[8]) != t_INT) mat[2]=(long)MC; /* not a small nf */  /**                                                                **/
   mat[3]=(long)mulmat_real(MC,M);  /********************************************************************/
   /* remove duplicate polynomials in y, updating a (same indexes), in place */
   static void
   remove_duplicates(GEN y, GEN a)
   {
     long k, i, l = lg(y);
     gpmem_t av = avma;
     GEN z;
   
     if (l < 2) return;
     z = new_chunk(3);
     z[1] = (long)y;
     z[2] = (long)a; (void)sort_factor(z, gcmp);
     for  (k=1, i=2; i<l; i++)
       if (!gegal((GEN)y[i], (GEN)y[k]))
       {
         k++;
         a[k] = a[i];
         y[k] = y[i];
       }
     l = k+1; setlg(a,l); setlg(y,l);
     avma = av;
   }
   
   /* if CHECK != NULL, return the first polynomial pol found such that
    * CHECK->f(CHECK->data, pol) != 0; return NULL if there are no such pol */
   static GEN
   _polred(GEN x, GEN a, GEN *pta, FP_chk_fun *CHECK)
   {
     long i, v = varn(x), l = lg(a);
     GEN ch, d, y = cgetg(l,t_VEC);
   
     for (i=1; i<l; i++)
     {
       if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); }
       ch = QX_caract(x, (GEN)a[i], v);
       if (CHECK)
       {
         ch = CHECK->f(CHECK->data, ch);
         if (!ch) continue;
         return ch;
       }
       d = modulargcd(derivpol(ch), ch);
       if (degpol(d)) ch = gdivexact(ch,d);
   
       if (canon_pol(ch) < 0 && pta) a[i] = (long)gneg_i((GEN)a[i]);
       if (DEBUGLEVEL > 3) outerr(ch);
       y[i] = (long)ch;
     }
     if (CHECK) return NULL; /* no suitable polynomial found */
     remove_duplicates(y,a);
     if (pta) *pta = a;
     return y;
   }
   
   static GEN
   allpolred(GEN x, long flag, GEN fa, GEN *pta, FP_chk_fun *CHECK)
   {
     GEN ro = NULL;
     nfbasic_t T; nfbasic_init(x, flag, fa, &T);
     if (T.lead) err(impl,"polred for non-monic polynomial");
     set_LLL_basis(&T, &ro);
     return _polred(T.x, T.bas, pta, CHECK);
   }
   
   /* FIXME: backward compatibility */
   #define red_PARTIAL 1
   #define red_ORIG    2
   
   GEN
   polred0(GEN x, long flag, GEN fa)
   {
     gpmem_t av = avma;
     GEN y, a;
     int fl = 0;
   
     if (fa && gcmp0(fa)) fa = NULL; /* compatibility */
     if (flag & red_PARTIAL) fl |= nf_PARTIALFACT;
     if (flag & red_ORIG)    fl |= nf_ORIG;
     y = allpolred(x, fl, fa, &a, NULL);
     if (fl & nf_ORIG) {
       GEN z = cgetg(3,t_MAT);
       z[1] = (long)a;
       z[2] = (long)y; y = z;
     }
   return gerepilecopy(av, y);    return gerepilecopy(av, y);
 }  }
   
   GEN
   ordred(GEN x)
   {
     gpmem_t av = avma;
     GEN y;
   
     if (typ(x) != t_POL) err(typeer,"ordred");
     if (!gcmp1(leading_term(x))) err(impl,"ordred");
     if (!signe(x)) return gcopy(x);
     y = cgetg(3,t_VEC);
     y[1] = (long)x;
     y[2] = (long)idmat(degpol(x));
     return gerepileupto(av, polred(y));
   }
   
   GEN
   polredfirstpol(GEN x, long flag, FP_chk_fun *CHECK)
   { return allpolred(x,flag,NULL,NULL,CHECK); }
   GEN
   polred(GEN x)
   { return polred0(x, 0, NULL); }
   GEN
   smallpolred(GEN x)
   { return polred0(x, red_PARTIAL, NULL); }
   GEN
   factoredpolred(GEN x, GEN fa)
   { return polred0(x, 0, fa); }
   GEN
   polred2(GEN x)
   { return polred0(x, red_ORIG, NULL); }
   GEN
   smallpolred2(GEN x)
   { return polred0(x, red_PARTIAL|red_ORIG, NULL); }
   GEN
   factoredpolred2(GEN x, GEN fa)
   { return polred0(x, red_PARTIAL, fa); }
   
   /********************************************************************/
   /**                                                                **/
   /**                           POLREDABS                            **/
   /**                                                                **/
   /********************************************************************/
   
   GEN
   T2_from_embed_norm(GEN x, long r1)
   {
     GEN p = sum(x, 1, r1);
     GEN q = sum(x, r1+1, lg(x)-1);
     if (q != gzero) p = gadd(p, gmul2n(q,1));
     return p;
   }
   
   GEN
   T2_from_embed(GEN x, long r1)
   {
     return T2_from_embed_norm(gnorm(x), r1);
   }
   
   typedef struct {
     long r1, v, prec;
     GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
     GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
     GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
     GEN bound; /* T2 norm of the polynomial defining nf */
   } CG_data;
   
   /* characteristic pol of x */
   static GEN
   get_polchar(CG_data *d, GEN x)
   {
     GEN g = gmul(d->ZKembed,x);
     long e;
     g = grndtoi(roots_to_pol_r1r2(g, d->r1, d->v), &e);
     if (e > -5) err(precer, "get_polchar");
     return g;
   }
   
   /* return a defining polynomial for Q(x) */
   static GEN
   get_polmin(CG_data *d, GEN x)
   {
     GEN g = get_polchar(d,x);
     GEN h = modulargcd(g,derivpol(g));
     if (degpol(h)) g = gdivexact(g,h);
     return g;
   }
   
   /* does x generate the correct field ? */
   static GEN
   chk_gen(void *data, GEN x)
   {
     gpmem_t av = avma;
     GEN g = get_polchar((CG_data*)data,x);
     GEN h = modulargcd(g,derivpol(g));
     if (degpol(h)) { avma = av; return NULL; }
     if (DEBUGLEVEL>3) fprintferr("  generator: %Z\n",g);
     return g;
   }
   
   /* mat = base change matrix, gram = LLL-reduced T2 matrix */
   static GEN
   chk_gen_init(FP_chk_fun *chk, GEN gram, GEN mat)
   {
     CG_data *d = (CG_data*)chk->data;
     GEN P,bound,prev,x,B;
     long l = lg(gram), N = l-1,i,prec,prec2;
     int skipfirst = 0;
   
     d->u = mat;
     d->ZKembed = gmul(d->M,   mat);
   
     bound = d->bound;
     prev = NULL;
     x = zerocol(N);
     for (i=1; i<l; i++)
     {
       B = gcoeff(gram,i,i);
       if (gcmp(B,bound) >= 0) continue; /* don't assume increasing norms */
   
       x[i] = un; P = get_polmin(d,x);
       x[i] = zero;
       if (degpol(P) == N)
       {
         if (gcmp(B,bound) < 0) bound = B ;
         continue;
       }
       if (DEBUGLEVEL>2) fprintferr("chk_gen_init: subfield %Z\n",P);
       if (skipfirst != i-1) continue;
   
       /* Q(w[1],...,w[i-1]) = Q[X]/(prev) is a strict subfield of nf */
       if (prev && degpol(prev) < N && !gegal(prev,P))
       {
         if (degpol(prev) * degpol(P) > 64) continue; /* too expensive */
         P = compositum(prev,P);
         P = (GEN)P[lg(P)-1];
         if (degpol(P) == N) continue;
         if (DEBUGLEVEL>2 && degpol(P) != degpol(prev))
           fprintferr("chk_gen_init: subfield %Z\n",P);
       }
       skipfirst++; prev = P;
     }
     /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
     chk->skipfirst = skipfirst;
     if (DEBUGLEVEL>2) fprintferr("chk_gen_init: skipfirst = %ld\n",skipfirst);
   
     /* should be gexpo( max_k C^n_k (bound/k)^(k/2) ) */
     prec2 = (1 + (((gexpo(bound)*N)/2) >> TWOPOTBITS_IN_LONG));
     if (prec2 < 0) prec2 = 0;
     prec = 3 + prec2;
     if (DEBUGLEVEL)
       fprintferr("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
     if (prec > d->prec) err(bugparier, "precision problem in polredabs");
     if (prec < d->prec) d->ZKembed = gprec_w(d->ZKembed, prec);
     return bound;
   }
   
   /* store phi(beta mod z). Suitable for gerepileupto */
   static GEN
   storeeval(GEN beta, GEN z, GEN lead)
   {
     GEN y = cgetg(3,t_VEC);
     z = gcopy(z);
     beta = lead? gdiv(beta, lead): gcopy(beta);
     y[1] = (long)z;
     y[2] = (long)to_polmod(beta,z);
     return y;
   }
   
   static GEN
   storeraw(GEN beta, GEN z)
   {
     GEN y = cgetg(3,t_VEC);
     y[1] = lcopy(z);
     y[2] = lcopy(beta); return y;
   }
   
   static void
   findmindisc(GEN *py, GEN *pa)
   {
     GEN v, dmin, z, b, discs, y = *py, a = *pa;
     long i,k, l = lg(y);
   
     if (l == 2) { *py = (GEN)y[1]; *pa = (GEN)a[1]; return; }
   
     discs = cgetg(l,t_VEC);
     for (i=1; i<l; i++) discs[i] = labsi(ZX_disc((GEN)y[i]));
     v = sindexsort(discs);
     k = v[1];
     dmin = (GEN)discs[k];
     z    = (GEN)y[k];
     b = (GEN)a[k];
     for (i=2; i<l; i++)
     {
       k = v[i];
       if (!egalii((GEN)discs[k], dmin)) break;
       if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; b = (GEN)a[k]; }
     }
     *py = z; *pa = b;
   }
   
   static GEN
   storepol(GEN x, GEN z, GEN a, GEN lead, long flag)
   {
     GEN y, b;
     if (flag & nf_RAW)
       y = storeraw(a, z);
     else if (flag & nf_ORIG)
     {
       b = modreverse_i(a, x);
       y = storeeval(b,z, lead);
     }
     else
       y = gcopy(z);
     return y;
   }
   
   static GEN
   storeallpol(GEN x, GEN z, GEN a, GEN lead, long flag)
   {
     GEN y, b;
   
     if (flag & nf_RAW)
     {
       long i, c = lg(z);
       y = cgetg(c,t_VEC);
       for (i=1; i<c; i++) y[i] = (long)storeraw((GEN)a[i], (GEN)z[i]);
     }
     else if (flag & nf_ORIG)
     {
       long i, c = lg(z);
       b = cgetg(c, t_VEC);
       for (i=1; i<c; i++)
         b[i] = (long)modreverse_i((GEN)a[i], x);
   
       y = cgetg(c,t_VEC);
       for (i=1; i<c; i++) y[i] = (long)storeeval((GEN)b[i], (GEN)z[i], lead);
     }
     else
       y = gcopy(z);
     return y;
   }
   
   GEN
   _polredabs(nfbasic_t *T, GEN *u)
   {
     long i, prec, e, n = degpol(T->x);
     GEN v, ro = NULL;
     FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, 0 };
     nffp_t F;
     CG_data d; chk.data = (void*)&d;
   
     set_LLL_basis(T, &ro);
     if (DEBUGLEVEL) msgtimer("LLL basis");
   
     /* || polchar ||_oo < 2^e */
     e = n * (gexpo(gmulsg(n, cauchy_bound(T->x))) + 1);
     prec = DEFAULTPREC + (e >> TWOPOTBITS_IN_LONG);
   
     get_nf_fp_compo(T, &F, ro, prec);
   
     d.v = varn(T->x);
     d.r1= T->r1;
     d.bound = gmul(T2_from_embed(F.ro, d.r1), dbltor(1.00000001));
     for (i=1; ; i++)
     {
       GEN R = R_from_QR(F.G, prec);
       d.prec = prec;
       d.M    = F.M;
       if (R)
       {
         v = fincke_pohst(_vec(R),NULL,5000,1, 0, &chk);
         if (v) break;
       }
       if (i==MAXITERPOL) err(accurer,"polredabs0");
       prec = (prec<<1)-2;
       get_nf_fp_compo(T, &F, NULL, prec);
       if (DEBUGLEVEL) err(warnprec,"polredabs0",prec);
     }
     *u = d.u; return v;
   }
   
   GEN
   polredabs0(GEN x, long flag)
   {
     long i, l, vx;
     gpmem_t av = avma;
     GEN y, a, u;
     nfbasic_t T;
   
     nfbasic_init(x, flag & nf_PARTIALFACT, NULL, &T);
     x = T.x; vx = varn(x);
   
     if (degpol(x) == 1)
     {
       u = NULL;
       y = _vec(polx[vx]);
       a = _vec(gsub((GEN)y[1], (GEN)x[2]));
     }
     else
     {
       GEN v = _polredabs(&T, &u);
       y = (GEN)v[1];
       a = (GEN)v[2];
     }
     l = lg(a);
     for (i=1; i<l; i++)
       if (canon_pol((GEN)y[i]) < 0) a[i] = (long)gneg_i((GEN)a[i]);
     remove_duplicates(y,a);
     l = lg(a);
     if (l == 1)
     {
       y = _vec(x);
       a = _vec(polx[vx]);
     }
     if (DEBUGLEVEL) fprintferr("%ld minimal vectors found.\n",l-1);
     if (flag & nf_ALL)
     {
       if (u) for (i=1; i<l; i++) a[i] = lmul(T.bas, gmul(u, (GEN)a[i]));
       y = storeallpol(x,y,a,T.lead,flag);
     }
     else
     {
       findmindisc(&y, &a);
       if (u) a = gmul(T.bas, gmul(u, a));
       y = storepol(x,y,a,T.lead,flag);
     }
     return gerepileupto(av, y);
   }
   
   GEN
   polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
   GEN
   polredabs(GEN x) { return polredabs0(x,0); }
   GEN
   polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
   
 static long  static long
 nf_pm1(GEN y)  nf_pm1(GEN y)
 {  {
Line 1340  nf_pm1(GEN y)
Line 1998  nf_pm1(GEN y)
   for (i=2; i<l; i++)    for (i=2; i<l; i++)
     if (signe(y[i])) return 0;      if (signe(y[i])) return 0;
   return signe(y[1]);    return signe(y[1]);
   
 }  }
   
 static GEN  static GEN
Line 1362  is_primitive_root(GEN nf, GEN fa, GEN x, long w)
Line 2019  is_primitive_root(GEN nf, GEN fa, GEN x, long w)
   return x;    return x;
 }  }
   
   /* only roots of 1 are +/- 1 */
   static GEN
   trivroots(GEN nf)
   {
     GEN y = cgetg(3, t_VEC);
     y[1] = deux;
     y[2] = (long)gscalcol_i(negi(gun), degpol(nf[1]));
     return y;
   }
   
 GEN  GEN
 rootsof1(GEN nf)  rootsof1(GEN nf)
 {  {
   ulong av;    gpmem_t av = avma;
   long N,k,i,ws,prec;    long N,k,i,ws,prec;
   GEN algun,p1,y,R1,d,list,w;    GEN p1,y,d,list,w;
   
   y=cgetg(3,t_VEC); av=avma; nf=checknf(nf);    nf = checknf(nf);
   R1=gmael(nf,2,1); algun=gmael(nf,8,1);    if ( nf_get_r1(nf) ) return trivroots(nf);
   if (signe(R1))  
   {    N = degpol(nf[1]); prec = nfgetprec(nf);
     y[1]=deux;  
     y[2]=lneg(algun); return y;  
   }  
   N=degpol(nf[1]); prec=gprecision((GEN)nf[6]);  
 #ifdef LONG_IS_32BIT  
   if (prec < 10) prec = 10;  
 #else  
   if (prec < 6) prec = 6;  
 #endif  
   for (i=1; ; i++)    for (i=1; ; i++)
   {    {
     p1 = fincke_pohst(nf,stoi(N),1000,1,prec,NULL);      GEN R = R_from_QR(gmael(nf,5,2), prec);
     if (p1) break;      if (R)
       {
         p1 = fincke_pohst(_vec(R),stoi(N),1000,1, 0, NULL);
         if (p1) break;
       }
     if (i == MAXITERPOL) err(accurer,"rootsof1");      if (i == MAXITERPOL) err(accurer,"rootsof1");
     prec=(prec<<1)-2;      prec = (prec<<1)-2;
     if (DEBUGLEVEL) err(warnprec,"rootsof1",prec);      if (DEBUGLEVEL) err(warnprec,"rootsof1",prec);
     nf=nfnewprec(nf,prec);      nf = nfnewprec(nf,prec);
   }    }
   if (itos(ground((GEN)p1[2])) != N) err(bugparier,"rootsof1 (bug1)");    if (itos(ground((GEN)p1[2])) != N) err(bugparier,"rootsof1 (bug1)");
   w=(GEN)p1[1]; ws = itos(w);    w = (GEN)p1[1]; ws = itos(w);
   if (ws == 2)    if (ws == 2) { avma = av; return trivroots(nf); }
   {  
     y[1]=deux; avma=av;  
     y[2]=lneg(algun); return y;  
   }  
   
   d = decomp(w); list = (GEN)p1[3]; k = lg(list);    d = decomp(w); list = (GEN)p1[3]; k = lg(list);
   for (i=1; i<k; i++)    for (i=1; i<k; i++)
   {    {
     p1 = (GEN)list[i];      p1 = (GEN)list[i];
     p1 = is_primitive_root(nf,d,p1,ws);      p1 = is_primitive_root(nf,d,p1,ws);
     if (p1)      if (p1) break;
     {  
       y[2]=lpilecopy(av,p1);  
       y[1]=lstoi(ws); return y;  
     }  
   }    }
   err(bugparier,"rootsof1");    if (!p1) err(bugparier,"rootsof1");
   return NULL; /* not reached */    y = cgetg(3, t_VEC);
     y[1] = lstoi(ws);
     y[2] = (long)p1; return gerepilecopy(av, y);
 }  }
   
 /*******************************************************************/  /*******************************************************************/
Line 1419  rootsof1(GEN nf)
Line 2075  rootsof1(GEN nf)
 /*                     DEDEKIND ZETA FUNCTION                      */  /*                     DEDEKIND ZETA FUNCTION                      */
 /*                                                                 */  /*                                                                 */
 /*******************************************************************/  /*******************************************************************/
   
 ulong smulss(ulong x, ulong y, ulong *rem);  
   
 static GEN  static GEN
 dirzetak0(GEN nf, long N0)  dirzetak0(GEN nf, long N0)
 {  {
   GEN vect,p1,pol,disc,c,c2;    GEN vect,p1,pol,disc,c,c2;
   long av=avma,i,j,k,limk,lx;    gpmem_t av=avma;
     long i,j,k,limk,lx;
   ulong q,p,rem;    ulong q,p,rem;
   byteptr d=diffptr;    byteptr d=diffptr;
   long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0};    long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
   
   pol=(GEN)nf[1]; disc=(GEN)nf[4];    pol=(GEN)nf[1]; disc=(GEN)nf[4];
   c  = (GEN) gpmalloc((N0+1)*sizeof(long));    c  = (GEN) gpmalloc((N0+1)*sizeof(long));
Line 1440  dirzetak0(GEN nf, long N0)
Line 2094  dirzetak0(GEN nf, long N0)
   
   while (court[2]<=N0)    while (court[2]<=N0)
   {    {
     court[2] += *d++; if (! *d) err(primer1);      NEXT_PRIME_VIADIFF_CHECK(court[2], d);
     if (smodis(disc,court[2])) /* court does not divide index */      if (smodis(disc,court[2])) /* court does not divide index */
       { vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); }        { vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); }
     else      else
Line 1494  initzeta(GEN pol, long prec)
Line 2148  initzeta(GEN pol, long prec)
   GEN nfz,nf,alpha,beta,mu,gr1,gr2,gru,p1,p2,cst,A0,c0,c1,c2,eps,coef;    GEN nfz,nf,alpha,beta,mu,gr1,gr2,gru,p1,p2,cst,A0,c0,c1,c2,eps,coef;
   GEN limx,bnf,resi,zet,C,coeflog,racpi,aij,tabj,colzero, *tabcstn, *tabcstni;    GEN limx,bnf,resi,zet,C,coeflog,racpi,aij,tabj,colzero, *tabcstn, *tabcstni;
   GEN c_even,ck_even,c_odd,ck_odd,serie_even,serie_odd,serie_exp,Pi;    GEN c_even,ck_even,c_odd,ck_odd,serie_even,serie_odd,serie_exp,Pi;
   long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n, av,av2,tetpil;    long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n;
   long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0};    gpmem_t av,av2,tetpil;
     long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
   stackzone *zone, *zone0, *zone1;    stackzone *zone, *zone0, *zone1;
   
   /*************** Calcul du residu et des constantes ***************/    /*************** Calcul du residu et des constantes ***************/
   eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5);    eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5);
   nfz=cgetg(10,t_VEC);    nfz=cgetg(10,t_VEC);
   bnf=buchinit(pol,p1,p1,prec+1); prec=(prec<<1)-1;    bnf = bnfinit0(pol,2,NULL,prec+1); prec=(prec<<1)-1;
   bnf = checkbnf_discard(bnf);    bnf = checkbnf_discard(bnf);
   Pi = mppi(prec); racpi=gsqrt(Pi,prec);    Pi = mppi(prec); racpi=gsqrt(Pi,prec);
   
Line 1540  initzeta(GEN pol, long prec)
Line 2195  initzeta(GEN pol, long prec)
   if (is_bigint(p1) || N0 > 10000000)    if (is_bigint(p1) || N0 > 10000000)
     err(talker,"discriminant too large for initzeta, sorry");      err(talker,"discriminant too large for initzeta, sorry");
   if (DEBUGLEVEL>=2)    if (DEBUGLEVEL>=2)
     { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); timer2(); }      { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); (void)timer2(); }
   
   /* Calcul de imax */    /* Calcul de imax */
   
Line 1563  initzeta(GEN pol, long prec)
Line 2218  initzeta(GEN pol, long prec)
   zone  = switch_stack(NULL,i + 2*(N0+1) + 6*prec);    zone  = switch_stack(NULL,i + 2*(N0+1) + 6*prec);
   zone1 = switch_stack(NULL,2*i);    zone1 = switch_stack(NULL,2*i);
   zone0 = switch_stack(NULL,2*i);    zone0 = switch_stack(NULL,2*i);
   switch_stack(zone,1);    (void)switch_stack(zone,1);
   tabcstn  = (GEN*) cgetg(N0+1,t_VEC);    tabcstn  = (GEN*) cgetg(N0+1,t_VEC);
   tabcstni = (GEN*) cgetg(N0+1,t_VEC);    tabcstni = (GEN*) cgetg(N0+1,t_VEC);
   p1 = ginv(cst);    p1 = ginv(cst);
   for (i=1; i<=N0; i++) { tabcstni[i] = gun; tabcstn[i] = mulsr(i,p1); }    for (i=1; i<=N0; i++) { tabcstni[i] = gun; tabcstn[i] = mulsr(i,p1); }
   switch_stack(zone,0);    (void)switch_stack(zone,0);
   
   /********** Calcul des coefficients a(i,j) independants de s **********/    /********** Calcul des coefficients a(i,j) independants de s **********/
   
Line 1595  initzeta(GEN pol, long prec)
Line 2250  initzeta(GEN pol, long prec)
     if (k&1) ck_odd[k]=lneg((GEN)ck_odd[k]);      if (k&1) ck_odd[k]=lneg((GEN)ck_odd[k]);
     ck_odd[k]=ladd(gru,(GEN)ck_odd[k]);      ck_odd[k]=ladd(gru,(GEN)ck_odd[k]);
   }    }
   ck_odd[1]=lsub((GEN)ck_odd[1],gmul(gr1,glog(gdeux,prec)));    ck_odd[1]=lsub((GEN)ck_odd[1],gmul(gr1,mplog2(prec)));
   serie_even =cgetg(ru+3,t_SER); serie_odd=cgetg(r2+3,t_SER);    serie_even =cgetg(ru+3,t_SER); serie_odd=cgetg(r2+3,t_SER);
   serie_even[1] = serie_odd[1] = evalsigne(1)+evalvalp(1);    serie_even[1] = serie_odd[1] = evalsigne(1)+evalvalp(1);
   i=0;    i=0;
Line 1656  initzeta(GEN pol, long prec)
Line 2311  initzeta(GEN pol, long prec)
       t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1);        t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1);
       for (j=2; j<=ru; j++)        for (j=2; j<=ru; j++)
       {        {
         long av2 = avma; p2 = gmul((GEN)t[j],p1);          gpmem_t av2 = avma;
           p2 = gmul((GEN)t[j], p1);
         t[j+1] = (long)gerepileuptoleaf(av2, divrs(p2,j));          t[j+1] = (long)gerepileuptoleaf(av2, divrs(p2,j));
       }        }
       /* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */        /* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */
Line 1705  initzeta(GEN pol, long prec)
Line 2361  initzeta(GEN pol, long prec)
     }      }
     /* use a parallel stack */      /* use a parallel stack */
     z = i&1? zone1: zone0;      z = i&1? zone1: zone0;
     switch_stack(z, 1);      (void)switch_stack(z, 1);
     for (n=1; n<=N0; n++)      for (n=1; n<=N0; n++)
       if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]);        if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]);
     /* come back */      /* come back */
     switch_stack(z, 0);      (void)switch_stack(z, 0);
   }    }
   nfz[4] = (long) C;    nfz[4] = (long) C;
   if (DEBUGLEVEL>=2) msgtimer("Cik");    if (DEBUGLEVEL>=2) msgtimer("Cik");
Line 1724  gzetakall(GEN nfz, GEN s, long flag, long prec2)
Line 2380  gzetakall(GEN nfz, GEN s, long flag, long prec2)
   GEN resi,C,cst,cstlog,coeflog,cs,coef;    GEN resi,C,cst,cstlog,coeflog,cs,coef;
   GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2;    GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2;
   GEN p1,unmoins,gexpro,gar,val,valm,valk,valkm;    GEN p1,unmoins,gexpro,gar,val,valm,valk,valkm;
   long ts = typ(s), r1,r2,ru,imax,i,j,k,N0,sl,prec,bigprec, av = avma;    long ts = typ(s), r1,r2,ru,imax,i,j,k,N0,sl,prec,bigprec;
     gpmem_t av = avma;
   
   if (typ(nfz)!=t_VEC || lg(nfz)!=10 || typ(nfz[1]) != t_VEC)    if (typ(nfz)!=t_VEC || lg(nfz)!=10 || typ(nfz[1]) != t_VEC)
     err(talker,"not a zeta number field in zetakall");      err(talker,"not a zeta number field in zetakall");
Line 1824  gzetakall(GEN nfz, GEN s, long flag, long prec2)
Line 2481  gzetakall(GEN nfz, GEN s, long flag, long prec2)
         for (k=1; k<=ru; k++)          for (k=1; k<=ru; k++)
         {          {
             p1 = gcoeff(C,i,k);              p1 = gcoeff(C,i,k);
             lambd = gsub(lambd,gdiv(p1,valk )); valk  = mulii(val,valk);              lambd = gsub(lambd,gdiv(p1,valk )); valk  = mulii(val,valk);
             lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);              lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);
         }          }
         val  = addis(val,1);          val  = addis(val,1);
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>