===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/base1.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/base1.c	2001/10/02 11:17:01	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/base1.c	2002/09/11 07:26:48	1.2
@@ -1,4 +1,4 @@
-/* $Id: base1.c,v 1.1 2001/10/02 11:17:01 noro Exp $
+/* $Id: base1.c,v 1.2 2002/09/11 07:26:48 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -20,35 +20,66 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /**************************************************************/
 #include "pari.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
 checkrnf(GEN rnf)
 {
-  if (typ(rnf)!=t_VEC) err(idealer1);
-  if (lg(rnf)!=12) err(idealer1);
+  if (typ(rnf)!=t_VEC || lg(rnf)!=12) err(idealer1);
 }
 
 GEN
-checkbnf(GEN bnf)
+_checkbnf(GEN bnf)
 {
-  if (typ(bnf)!=t_VEC) err(idealer1);
-  switch (lg(bnf))
+  if (typ(bnf) == t_VEC)
+    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;
-    case 10:
-      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]);
+    if (_checknf(x)) err(talker,"please apply bnfinit first");
+    err(idealer1);
   }
-  err(idealer1);
-  return NULL; /* not reached */
+  return bnf;
 }
 
 GEN
@@ -61,19 +92,15 @@ checkbnf_discard(GEN bnf)
 }
 
 GEN
-checknf(GEN nf)
+checknf(GEN x)
 {
-  if (typ(nf)==t_POL) err(talker,"please apply nfinit first");
-  if (typ(nf)!=t_VEC) err(idealer1);
-  switch(lg(nf))
+  GEN nf = _checknf(x);
+  if (!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]);
+    if (typ(x)==t_POL) err(talker,"please apply nfinit first");
+    err(idealer1);
   }
-  err(idealer1);
-  return NULL; /* not reached */
+  return nf;
 }
 
 void
@@ -126,13 +153,6 @@ checkprimeid(GEN id)
 }
 
 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)
 {
   long k = lgef(x)-1;
@@ -248,7 +268,8 @@ transroot(GEN x, int i, int j)
 GEN
 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);
 
   if (typ(x)!=t_POL) err(notpoler,"tschirnhaus");
@@ -260,12 +281,12 @@ tschirnhaus(GEN x)
     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[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 */
   if (DEBUGLEVEL>1)
     fprintferr("Tschirnhaus transform. New pol: %Z",u);
-  avma=a; return gerepileupto(av,u);
+  avma=av2; return gerepileupto(av,u);
 }
 #undef randshift
 
@@ -284,7 +305,8 @@ gpolcomp(GEN p1, GEN p2)
   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
 primitive_pol_to_monic(GEN pol, GEN *ptlead)
 {
@@ -342,14 +364,32 @@ get_F4(GEN x)
   return p1;
 }
 
-GEN galoisbig(GEN x, long prec);
-GEN roots_to_pol(GEN a, long v);
+static GEN
+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
 galois(GEN x, long prec)
 {
-  long av=avma,av1,i,j,k,n,f,l,l2,e,e1,pr,ind;
-  GEN x1,p1,p2,p3,p4,p5,p6,w,y,z,ee;
+  gpmem_t av = avma, av1;
+  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 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,
@@ -357,35 +397,29 @@ galois(GEN x, long prec)
   if (typ(x)!=t_POL) err(notpoler,"galois");
   n=degpol(x); if (n<=0) err(constpoler,"galois");
   if (n>11) err(impl,"galois of degree higher than 11");
-  x = gdiv(x,content(x));
+  x = primpart(x);
   check_pol_int(x, "galois");
   if (gisirreducible(x) != gun)
     err(impl,"galois of reducible polynomial");
 
   if (n<4)
   {
-    if (n<3)
-    {
-      avma=av; y=cgetg(4,t_VEC);
-      y[1] = (n==1)? un: deux;
-      y[2]=lnegi(gun);
-    }
-    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;
+    if (n == 1) { avma = av; return _res(1, 1,1); }
+    if (n == 2) { avma = av; return _res(2,-1,1); }
+    /* n = 3 */
+    f = carreparfait(ZX_disc(x));
+    avma = av;
+    return f? _res(3,1,1): _res(6,-1,2);
   }
   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(;;)
   {
+    GEN cb = cauchy_bound(x);
     switch(n)
     {
       case 4: z = cgetg(7,t_VEC);
+        prec = DEFAULTPREC + (long)((gexpo(cb)*18. / BITS_IN_LONG));
         for(;;)
 	{
 	  p1=roots(x,prec);
@@ -395,35 +429,29 @@ galois(GEN x, long prec)
           z[4] = (long)get_F4(transroot(p1,1,4));
           z[5] = (long)get_F4(transroot(p1,2,3));
           z[6] = (long)get_F4(transroot(p1,3,4));
-          p4 = roots_to_pol(z,0);
-          p5=grndtoi(greal(p4),&e);
-          e1=gexpo(gimag(p4)); if (e1>e) e=e1;
+          p5 = roots_to_ZX(z, &e);
           if (e <= -10) break;
 	  prec = (prec<<1)-2;
 	}
-	p6 = modulargcd(p5, derivpol(p5));
-	if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
+        if (!ZX_is_squarefree(p5)) goto tchi;
 	p2 = (GEN)factor(p5)[1];
 	switch(lg(p2)-1)
 	{
-	  case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);
-	    y[3]=un;
-	    if (f) { y[2]=un; y[1]=lstoi(12); return y; }
-	    y[2]=lnegi(gun); y[1]=lstoi(24); return y;
+	  case 1: f = carreparfait(ZX_disc(x)); avma = av;
+            return f? _res(12,1,4): _res(24,-1,5);
 
-	  case 2: avma=av; y=cgetg(4,t_VEC);
-	    y[3]=un; y[2]=lnegi(gun); y[1]=lstoi(8); return y;
+	  case 2: avma = av; return _res(8,-1,3);
 	
-	  case 3: avma=av; y=cgetg(4,t_VEC);
-	    y[1]=lstoi(4); y[3]=un;
-	    y[2] = (lgef(p2[1])==5)? un: lnegi(gun);
-	    return y;
+	  case 3: avma = av;
+	    return (degpol(p2[1])==2)? _res(4,1,2): _res(4,-1,1);
 
 	  default: err(bugparier,"galois (bug1)");
 	}
 
       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(;;)
@@ -443,27 +471,20 @@ galois(GEN x, long prec)
               e1 = gexpo(gimag(p3)); if (e1>e) e=e1;
               ee[l]=e; z[l] = (long)p3;
 	    }
-            p4 = roots_to_pol(z,0);
-	    p5=grndtoi(greal(p4),&e);
-            e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
+            p5 = roots_to_ZX(z, &e);
             if (e <= -10) break;
 	    prec = (prec<<1)-2;
 	  }
-	  p6 = modulargcd(p5,derivpol(p5));
-	  if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
+          if (!ZX_is_squarefree(p5)) goto tchi;
 	  p3=(GEN)factor(p5)[1];
-	  f=carreparfait(discsr(x));
+	  f=carreparfait(ZX_disc(x));
 	  if (lg(p3)-1==1)
 	  {
-	    avma=av; y=cgetg(4,t_VEC); y[3]=un;
-	    if (f) { y[2]=un; y[1]=lstoi(60); return y; }
-	    else { y[2]=lneg(gun); y[1]=lstoi(120); return y; }
+	    avma = av;
+            return f? _res(60,1,4): _res(120,-1,5);
 	  }
-	  if (!f)
-	  {
-	    avma=av; y=cgetg(4,t_VEC);
-	    y[3]=un; y[2]=lneg(gun); y[1]=lstoi(20); return y;
-	  }
+	  if (!f) { avma = av; return _res(20,-1,3); }
+
           pr = - (bit_accuracy(prec) >> 1);
           for (l=1; l<=6; l++)
 	    if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break;
@@ -472,23 +493,23 @@ galois(GEN x, long prec)
 	  p3=gzero;
 	  for (i=1; i<=5; i++)
 	  {
-	    j=(i%5)+1;
-	    p3=gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]),
-			    gsub((GEN)p2[j],(GEN)p2[i])));
+	    j = i+1; if (j>5) j -= 5;
+	    p3 = gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]),
+		              gsub((GEN)p2[j],(GEN)p2[i])));
 	  }
 	  p5=gsqr(p3); p4=grndtoi(greal(p5),&e);
           e1 = gexpo(gimag(p5)); if (e1>e) e=e1;
 	  if (e <= -10)
 	  {
 	    if (gcmp0(p4)) goto tchi;
-	    f=carreparfait(p4); avma=av; y=cgetg(4,t_VEC);
-	    y[3]=y[2]=un; y[1]=lstoi(f?5:10);
-	    return y;
+	    f = carreparfait(p4); avma = av;
+            return f? _res(5,1,1): _res(10,1,2);
 	  }
 	  prec=(prec<<1)-2;
 	}
 
       case 6: z = cgetg(7, t_VEC);
+        prec = DEFAULTPREC + (long)((gexpo(cb)*42. / BITS_IN_LONG));
         for(;;)
 	{
           for(;;)
@@ -502,18 +523,16 @@ galois(GEN x, long prec)
 	      {
 		p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]),
 		        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;
 	    }
-            p4=roots_to_pol(z,0);
-	    p5=grndtoi(greal(p4),&e);
-            e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
+            p5 = roots_to_ZX(z, &e);
             if (e <= -10) break;
 	    prec=(prec<<1)-2;
 	  }
-	  p6 = modulargcd(p5,derivpol(p5));
-	  if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
+	  if (!ZX_is_squarefree(p5)) goto tchi;
 	  p2=(GEN)factor(p5)[1];
 	  switch(lg(p2)-1)
 	  {
@@ -530,127 +549,93 @@ galois(GEN x, long prec)
 		          gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6]));
 		  z[++ind] = (long)p3;
 		}
-              p4 = roots_to_pol(z,0);
-	      p5=grndtoi(greal(p4),&e);
-              e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
+              p5 = roots_to_ZX(z, &e);
 	      if (e <= -10)
 	      {
-		p6 = modulargcd(p5,derivpol(p5));
-		if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
-		p2=(GEN)factor(p5)[1];
-		f=carreparfait(discsr(x));
-		avma=av; y=cgetg(4,t_VEC); y[3]=un;
+                if (!ZX_is_squarefree(p5)) goto tchi;
+		p2 = (GEN)factor(p5)[1];
+		f = carreparfait(ZX_disc(x));
+		avma = av;
 		if (lg(p2)-1==1)
-		{
-		  if (f) { y[2]=un; y[1]=lstoi(360); }
-		  else { y[2]=lnegi(gun); y[1]=lstoi(720); }
-		}
+                  return f? _res(360,1,15): _res(720,-1,16);
 		else
-		{
-		  if (f) { y[2]=un; y[1]=lstoi(36); }
-		  else { y[2]=lnegi(gun); y[1]=lstoi(72); }
-		}
-                return y;
+                  return f? _res(36,1,10): _res(72,-1,13);
 	      }
 	      prec=(prec<<1)-2; break;
 		
 	    case 2: l2=degpol(p2[1]); if (l2>3) l2=6-l2;
 	      switch(l2)
 	      {
-		case 1: f=carreparfait(discsr(x));
-		  avma=av; y=cgetg(4,t_VEC); y[3]=un;
-		  if (f) { y[2]=un; y[1]=lstoi(60); }
-		  else { y[2]=lneg(gun); y[1]=lstoi(120); }
-		  return y;
-		case 2: f=carreparfait(discsr(x));
-		  if (f)
-		  {
-		    avma=av; y=cgetg(4,t_VEC);
-		    y[3]=y[2]=un; y[1]=lstoi(24);
-		  }
-		  else
-		  {
-		    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 1: f = carreparfait(ZX_disc(x));
+		  avma = av;
+                  return f? _res(60,1,12): _res(120,-1,14);
+		case 2: f = carreparfait(ZX_disc(x));
+		  if (f) { avma = av; return _res(24,1,7); }
+                  p3 = (degpol(p2[1])==2)? (GEN)p2[2]: (GEN)p2[1];
+                  f = carreparfait(ZX_disc(p3));
+                  avma = av;
+                  return f? _res(24,-1,6): _res(48,-1,11);
+		case 3: f = carreparfait(ZX_disc((GEN)p2[1]))
+		         || carreparfait(ZX_disc((GEN)p2[2]));
+		  avma = av;
+                  return f? _res(18,-1,5): _res(36,-1,9);
 	      }
 	    case 3:
 	      for (l2=1; l2<=3; l2++)
-		if (lgef(p2[l2])>=6) p3=(GEN)p2[l2];
-	      if (lgef(p3)==6)
+		if (degpol(p2[l2]) >= 3) p3 = (GEN)p2[l2];
+	      if (degpol(p3) == 3)
 	      {
-		f=carreparfait(discsr(p3)); avma=av; y=cgetg(4,t_VEC);
-                y[2]=lneg(gun); y[1]=lstoi(f? 6: 12);
+		f = carreparfait(ZX_disc(p3)); avma = av;
+                return f? _res(6,-1,1): _res(12,-1,3);
 	      }
 	      else
 	      {
-		f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);
-		if (f) { y[2]=un; y[1]=lstoi(12); }
-		else { y[2]=lneg(gun); y[1]=lstoi(24); }
+		f = carreparfait(ZX_disc(x)); avma = av;
+                return f? _res(12,1,4): _res(24,-1,8);
 	      }
-              y[3]=un; return y;
-	    case 4: avma=av; y=cgetg(4,t_VEC);
-	      y[1]=lstoi(6); y[2]=lneg(gun); y[3]=deux; return y;
+	    case 4: avma = av; return _res(6,-1,2);
             default: err(bugparier,"galois (bug3)");
 	  }
 	}
 	
       case 7: z = cgetg(36,t_VEC);
+        prec = DEFAULTPREC + (long)((gexpo(cb)*7. / BITS_IN_LONG));
         for(;;)
 	{
-	  ind = 0; p1=roots(x,prec); p4=gun;
+	  ind = 0; p1=roots(x,prec);
 	  for (i=1; i<=5; i++)
 	    for (j=i+1; j<=6; j++)
             {
-              p6 = gadd((GEN)p1[i],(GEN)p1[j]);
-	      for (k=j+1; k<=7; k++)
-                z[++ind] = (long) gadd(p6,(GEN)p1[k]);
+              GEN t = gadd((GEN)p1[i],(GEN)p1[j]);
+	      for (k=j+1; k<=7; k++) z[++ind] = ladd(t, (GEN)p1[k]);
             }
-          p4 = roots_to_pol(z,0);
-          p5 = grndtoi(greal(p4),&e);
-          e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
+          p5 = roots_to_ZX(z, &e);
 	  if (e <= -10) break;
           prec = (prec<<1)-2;
 	}
-	p6 = modulargcd(p5,derivpol(p5));
-	if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
+        if (!ZX_is_squarefree(p5)) goto tchi;
 	p2=(GEN)factor(p5)[1];
 	switch(lg(p2)-1)
 	{
-	  case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC); y[3]=un;
-	    if (f) { y[2]=un; y[1]=lstoi(2520); }
-	    else { y[2]=lneg(gun); y[1]=lstoi(5040); }
-	    return y;
-	  case 2: f=degpol(p2[1]); avma=av; y=cgetg(4,t_VEC); y[3]=un;
-	    if (f==7 || f==28) { y[2]=un; y[1]=lstoi(168); }
-	    else { y[2]=lneg(gun); y[1]=lstoi(42); }
-	    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;
+	  case 1: f = carreparfait(ZX_disc(x)); avma = av;
+            return f? _res(2520,1,6): _res(5040,-1,7);
+	  case 2: f = degpol(p2[1]); avma = av;
+	    return (f==7 || f==28)? _res(168,1,5): _res(42,-1,4);
+	  case 3: avma = av; return _res(21,1,3);
+	  case 4: avma = av; return _res(14,-1,2);
+	  case 5: avma = av; return _res(7,1,1);
           default: err(talker,"galois (bug2)");
 	}
     }
-    tchi: avma=av1; x=tschirnhaus(x1);
+    tchi: avma = av1; x = tschirnhaus(x1);
   }
 }
 
 GEN
 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;
 
   nf=checknf(nf); pol=(GEN)nf[1];
@@ -713,6 +698,16 @@ GEN pol_to_monic(GEN pol, GEN *lead);
 int cmp_pol(GEN x, GEN y);
 
 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)
 {
   if (typ(x) == t_POL) { *nf = NULL; return x; }
@@ -723,7 +718,7 @@ get_nfpol(GEN x, GEN *nf)
 static GEN
 nfiso0(GEN a, GEN b, long fliso)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long n,m,i,vb,lx;
   GEN nfa,nfb,p1,y,la,lb;
 
@@ -751,8 +746,8 @@ nfiso0(GEN a, GEN b, long fliso)
   }
   else
   {
-    GEN da = nfa? (GEN)nfa[3]: discsr(a);
-    GEN db = nfb? (GEN)nfb[3]: discsr(b);
+    GEN da = nfa? (GEN)nfa[3]: ZX_disc(a);
+    GEN db = nfb? (GEN)nfb[3]: ZX_disc(b);
     if (fliso)
     {
       p1=gdiv(da,db);
@@ -787,7 +782,7 @@ nfiso0(GEN a, GEN b, long fliso)
     }
     y = gen_sort(y, 0, cmp_pol);
     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; }
   for (i=1; i<lx; i++)
@@ -801,28 +796,28 @@ nfiso0(GEN a, GEN b, long fliso)
 }
 
 GEN
-nfisisom(GEN a, GEN b)
-{
-  return nfiso0(a,b,1);
-}
+nfisisom(GEN a, GEN b) { return nfiso0(a,b,1); }
 
 GEN
-nfisincl(GEN a, GEN b)
-{
-  return nfiso0(a,b,0);
-}
+nfisincl(GEN a, GEN b) { return nfiso0(a,b,0); }
 
 /*************************************************************************/
 /**									**/
 /**			       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) */
 GEN
 quicktrace(GEN x, GEN sym)
@@ -843,7 +838,7 @@ quicktrace(GEN x, GEN sym)
 static GEN
 trace_col(GEN x, GEN T)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN t = gzero;
   long i, l = lg(x);
 
@@ -853,92 +848,6 @@ trace_col(GEN x, GEN 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
  * as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing
  * by the content), and to to leading coeff otherwise.
@@ -948,84 +857,25 @@ GEN
 pol_to_monic(GEN pol, GEN *lead)
 {
   long n = lgef(pol)-1;
-  GEN p1;
 
   if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; }
-
-  p1=content(pol); if (!gcmp1(p1)) pol = gdiv(pol,p1);
-  return primitive_pol_to_monic(pol,lead);
+  return primitive_pol_to_monic(primpart(pol), lead);
 }
 
-/* NB: TI integral, det(TI) = d_K, ideal/dideal = codifferent */
-GEN
-make_MDI(GEN nf, GEN TI, GEN *ideal, GEN *dideal)
+/* Let T = (Tr(wi wj)), then T^(-1) Z^n = codifferent = coD/dcoD
+ * NB: det(T) = dK, TI = dK T^(-1) integral */
+static GEN
+make_MDI(GEN nf, GEN TI, GEN *coD, GEN *dcoD)
 {
-  GEN c = content(TI);
-  GEN p1;
+  GEN c, MDI, dK = (GEN)nf[3];
 
-  *dideal = divii((GEN)nf[3], c);
-  *ideal = p1 = hnfmodid(gdiv(TI,c), *dideal);
-  return gmul(c, ideal_two_elt(nf, p1));
+  TI = Q_primitive_part(TI, &c);
+  *dcoD = c? diviiexact(dK, c): dK;
+  *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] */
 long
 nf_get_r1(GEN nf)
@@ -1043,14 +893,22 @@ nf_get_r2(GEN nf)
     err(talker,"false nf in nf_get_r2");
   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
-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;
-  GEN tr,T,T1,sym;
   T = cgetg(n+1,t_MAT);
   T1 = cgetg(n+1,t_COL);
   sym = polsym(x, n-1);
@@ -1076,16 +934,17 @@ get_T(GEN mul, GEN x, GEN bas, GEN den)
 GEN
 get_bas_den(GEN bas)
 {
-  GEN z,d,den, dbas = dummycopy(bas);
-  long i, c = 0, l = lg(bas);
+  GEN z,b,d,den, dbas = dummycopy(bas);
+  long i, l = lg(bas);
+  int power = 1;
   den = cgetg(l,t_VEC);
   for (i=1; i<l; i++)
   {
-    d = denom(content((GEN)dbas[i]));
-    if (is_pm1(d)) d = NULL; else { dbas[i] = lmul((GEN)dbas[i],d); c++; }
-    den[i] = (long)d;
+    b = Q_remove_denom((GEN)bas[i], &d);
+    dbas[i]= (long)b;
+    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[1] = (long)dbas;
   z[2] = (long)den; return z;
@@ -1101,16 +960,16 @@ _mulii(GEN x, GEN y)
 }
 
 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);
   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 (i=1; i<=n; i++)
     for (j=i; j<=n; j++)
     {
-      ulong av = avma;
+      gpmem_t av = avma;
       z = gres(gmul((GEN)bas[j],(GEN)bas[i]), x);
       z = mulmat_pol(invbas, z); /* integral column */
       if (den)
@@ -1120,216 +979,1015 @@ get_mul_table(GEN x,GEN basden,GEN invbas,GEN *T)
       }
       mul[j+(i-1)*n] = mul[i+(j-1)*n] = lpileupto(av,z);
     }
-  if (T) *T = get_T(mul,x,bas,den);
   return mul;
 }
 
-/* fill mat = nf[5], as well as nf[8] and nf[9]
- * If (small) only compute a subset (use dummy 0s for the rest) */
+/* as get_Tr, mul_table not precomputed */
+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
-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];
-  GEN basden,mul,invbas,M,MC,T,MDI,D,TI,A,dA,mat;
-  long r1 = nf_get_r1(nf), n = lg(bas)-1;
+  F->x  = (GEN)nf[1];
+  F->ro = NULL;
+  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;
-  basden = get_bas_den(bas);
-  M = make_M(basden,ro);
-  MC = make_MC(r1,M);
-  mat[1]=(long)M;
-  mat[3]=(long)mulmat_real(MC,M);
-  if (small) { nf[8]=nf[9]=mat[2]=zero; return; }
+static void
+nffp_init(nffp_t *F, nfbasic_t *T, GEN ro, long prec)
+{
+  F->x  = T->x;
+  F->ro = ro;
+  F->r1 = T->r1;
+  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);
-  mul = get_mul_table(x,basden,invbas,&T);
+static void
+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");
-  
-  nf[8]=(long)invbas;
-  nf[9]=(long)mul;
 
-  TI = ZM_inv(T, dK); /* dK T^-1 */
-  MDI = make_MDI(nf,TI, &A, &dA);
-  mat[6]=(long)TI;
-  mat[7]=(long)MDI; /* needed in idealinv below */
-  if (is_pm1(index))
+  Tr = get_Tr((GEN)nf[9], x, F.basden);
+  TI = ZM_inv(Tr, T->dK); /* dK T^-1 */
+  mat[6] = (long)TI;
+  mat[7] = (long)make_MDI(nf,TI, &A, &dA); /* needed in idealinv below */
+  if (is_pm1(T->index))
     D = idealhermite_aux(nf, derivpol(x));
   else
     D = gmul(dA, idealinv(nf, A));
-  mat[2]=(long)MC;
-  mat[4]=(long)T;
-  mat[5]=(long)D;
+  mat[3] = zero; /* FIXME: was gram matrix of current mat[2]. Useless */
+  mat[4] = (long)Tr;
+  mat[5] = (long)D;
   if (DEBUGLEVEL) msgtimer("matrices");
+  return nf;
 }
 
-/* Initialize the number field defined by the polynomial x (in variable v)
- * flag & nf_REGULAR
- *    regular behaviour.
- * flag & nf_SMALL
- *    compute only nf[1] (pol), nf[2] (signature), nf[5][3] (T2) and
- *    nf[7] (integer basis), the other components are filled with gzero.
- * flag & nf_REDUCE
- *    try a polred first.
- * 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).
- */
+static GEN
+hnffromLLL(GEN nf)
+{
+  GEN d, x;
+  x = vecpol_to_mat((GEN)nf[7], degpol(nf[1]));
+  x = Q_remove_denom(x, &d);
+  if (!d) return x; /* power basis */
+  return gauss(hnfmodid(x, d), x);
+}
 
-/* 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
-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 av=avma,n,i,r1,r2,ru,PRECREG;
+  long tx = typ(x);
+  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();
-  if (typ(x)==t_POL)
+GEN
+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");
-    check_pol_int(x,"initalg");
-    if (gisirreducible(x) == gzero) err(redpoler,"nfinit");
-    if (!gcmp1((GEN)x[n+2]))
+    F.prec = prec; make_M_G(&F, 0); G = F.G;
+    if (u0) G = gmul(G, u0);
+    if (DEBUGLEVEL)
+      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 (!(flag & nf_SMALL))
-      {
-        if (!(flag & nf_REDUCE))
-          err(warner,"non-monic polynomial. Result of the form [nf,c]");
-        flag = flag | nf_REDUCE | nf_ORIG;
-      }
+      if (typ(u) == t_MAT) break;
+      u = (GEN)u[1];
+      if (u0) u0 = gerepileupto(av, gmul(u0,u));
+      else    u0 = gerepilecopy(av, u);
     }
-    bas = allbase4(x,0,&dK,NULL);
-    if (DEBUGLEVEL) msgtimer("round4");
-    dx = discsr(x); r1 = sturm(x);
+    if (i == MAXITERPOL) err(accurer,"red_T2");
+    prec = (prec<<1)-2 + (gexpo(u0) >> TWOPOTBITS_IN_LONG);
+    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
   {
-    i = lg(x);
-    if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)
-    { /* 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 */
+    u = lllint_marked(1, make_Tr(T->x, T->bas), 100, 1, &u,NULL,NULL);
+    if (!u) u = idmat(1);
   }
-  r2=(n-r1)>>1; ru=r1+r2;
-  PRECREG = prec + (expi(dK)>>(TWOPOTBITS_IN_LONG+1))
-                 + (long)((sqrt((double)n)+3) / sizeof(long) * 4);
-  rev = NULL;
-  if (flag & nf_REDUCE)
+  return u;
+}
+
+static void
+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);
-    if (DEBUGLEVEL) msgtimer("polred");
+    T->dxbest = dxn; T->xbest = xn; T->indbest = T->ind;
   }
-  if (!carrecomplet(divii(dx,dK),&index))
-    err(bugparier,"nfinit (incorrect discriminant)");
+  if (T->ind >= T->indmax) return T->xbest;
+  return NULL;
+}
 
-  ro=get_roots(x,r1,ru,PRECREG);
-  if (DEBUGLEVEL) msgtimer("roots");
+/* z in Z[X] with positive leading coeff. Set z := z(-X) or z(X) so that
+ * second coeff is > 0. In place. */
+static int
+canon_pol(GEN z)
+{
+  long i,s;
 
-  nf=cgetg(10,t_VEC);
-  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)
+  for (i=lgef(z)-2; i>=2; i-=2)
   {
-    if (!rev) err(talker,"bad flag in initalgall0");
-    res = cgetg(3,t_VEC);
-    res[1]=(long)nf; nf = res;
-    res[2]=lead? ldiv(rev,lead): (long)rev;
+    s = signe(z[i]);
+    if (s > 0)
+    {
+      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
-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
-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
+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)
 {
   switch(flag)
   {
     case 0:
-    case 1: return initalgall0(x,nf_REGULAR,prec);
-    case 2: return initalgall0(x,nf_REDUCE,prec);
-    case 3: return initalgall0(x,nf_REDUCE|nf_ORIG,prec);
-    case 4: return initalgall0(x,nf_REDUCE|nf_PARTIAL,prec);
-    case 5: return initalgall0(x,nf_REDUCE|nf_ORIG|nf_PARTIAL,prec);
-    case 6: return initalgall0(x,nf_SMALL,prec);
+    case 1: return _initalg(x,0,prec);
+    case 2: return _initalg(x,nf_RED,prec);
+    case 3: return _initalg(x,nf_RED|nf_ORIG,prec);
+    case 4: return _initalg(x,nf_PARTRED,prec);
+    case 5: return _initalg(x,nf_PARTRED|nf_ORIG,prec);
     default: err(flagerr,"nfinit");
   }
   return NULL; /* not reached */
 }
 
-GEN
-initalg(GEN x, long prec)
-{
-  return initalgall0(x,nf_REGULAR,prec);
-}
-
 /* assume x a bnr/bnf/nf */
 long
 nfgetprec(GEN x)
 {
   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
 nfnewprec(GEN nf, long prec)
 {
-  long av=avma,r1,r2,ru,n;
-  GEN y,pol,ro,basden,MC,mat,M;
+  gpmem_t av = avma;
+  GEN NF;
+  nffp_t F;
 
   if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec");
   if (lg(nf) == 11) return bnfnewprec(nf,prec);
   if (lg(nf) ==  7) return bnrnewprec(nf,prec);
   (void)checknf(nf);
-  y = dummycopy(nf);
-  pol=(GEN)nf[1]; n=degpol(pol);
-  r1 = nf_get_r1(nf);
-  r2 = (n - r1) >> 1; ru = r1+r2;
-  mat=dummycopy((GEN)nf[5]);
-  ro=get_roots(pol,r1,ru,prec);
-  y[5]=(long)mat;
-  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;
-  if (typ(nf[8]) != t_INT) mat[2]=(long)MC; /* not a small nf */
-  mat[3]=(long)mulmat_real(MC,M);
+  NF = dummycopy(nf);
+  NF[5] = (long)dummycopy((GEN)NF[5]);
+  remake_GM(NF, &F, prec);
+  NF[6]  = (long)F.ro;
+  mael(NF,5,1) = (long)F.M;
+  mael(NF,5,2) = (long)F.G;
+  return gerepilecopy(av, NF);
+}
+
+/********************************************************************/
+/**                                                                **/
+/**                           POLRED                               **/
+/**                                                                **/
+/********************************************************************/
+/* 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);
 }
 
+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
 nf_pm1(GEN y)
 {
@@ -1340,7 +1998,6 @@ nf_pm1(GEN y)
   for (i=2; i<l; i++)
     if (signe(y[i])) return 0;
   return signe(y[1]);
-
 }
 
 static GEN
@@ -1362,56 +2019,55 @@ is_primitive_root(GEN nf, GEN fa, GEN x, long w)
   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
 rootsof1(GEN nf)
 {
-  ulong av;
+  gpmem_t av = avma;
   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);
-  R1=gmael(nf,2,1); algun=gmael(nf,8,1);
-  if (signe(R1))
-  {
-    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
+  nf = checknf(nf);
+  if ( nf_get_r1(nf) ) return trivroots(nf);
+
+  N = degpol(nf[1]); prec = nfgetprec(nf);
   for (i=1; ; i++)
   {
-    p1 = fincke_pohst(nf,stoi(N),1000,1,prec,NULL);
-    if (p1) break;
+    GEN R = R_from_QR(gmael(nf,5,2), prec);
+    if (R)
+    {
+      p1 = fincke_pohst(_vec(R),stoi(N),1000,1, 0, NULL);
+      if (p1) break;
+    }
     if (i == MAXITERPOL) err(accurer,"rootsof1");
-    prec=(prec<<1)-2;
+    prec = (prec<<1)-2;
     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)");
-  w=(GEN)p1[1]; ws = itos(w);
-  if (ws == 2)
-  {
-    y[1]=deux; avma=av;
-    y[2]=lneg(algun); return y;
-  }
+  w = (GEN)p1[1]; ws = itos(w);
+  if (ws == 2) { avma = av; return trivroots(nf); }
 
   d = decomp(w); list = (GEN)p1[3]; k = lg(list);
   for (i=1; i<k; i++)
   {
     p1 = (GEN)list[i];
     p1 = is_primitive_root(nf,d,p1,ws);
-    if (p1)
-    {
-      y[2]=lpilecopy(av,p1);
-      y[1]=lstoi(ws); return y;
-    }
+    if (p1) break;
   }
-  err(bugparier,"rootsof1");
-  return NULL; /* not reached */
+  if (!p1) err(bugparier,"rootsof1");
+  y = cgetg(3, t_VEC);
+  y[1] = lstoi(ws);
+  y[2] = (long)p1; return gerepilecopy(av, y);
 }
 
 /*******************************************************************/
@@ -1419,17 +2075,15 @@ rootsof1(GEN nf)
 /*                     DEDEKIND ZETA FUNCTION                      */
 /*                                                                 */
 /*******************************************************************/
-
-ulong smulss(ulong x, ulong y, ulong *rem);
-
 static GEN
 dirzetak0(GEN nf, long N0)
 {
   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;
   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];
   c  = (GEN) gpmalloc((N0+1)*sizeof(long));
@@ -1440,7 +2094,7 @@ dirzetak0(GEN nf, long 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 */
       { vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); }
     else
@@ -1494,14 +2148,15 @@ initzeta(GEN pol, long prec)
   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 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 court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0};
+  long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n;
+  gpmem_t av,av2,tetpil;
+  long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
   stackzone *zone, *zone0, *zone1;
 
   /*************** Calcul du residu et des constantes ***************/
   eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5);
   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);
   Pi = mppi(prec); racpi=gsqrt(Pi,prec);
 
@@ -1540,7 +2195,7 @@ initzeta(GEN pol, long prec)
   if (is_bigint(p1) || N0 > 10000000)
     err(talker,"discriminant too large for initzeta, sorry");
   if (DEBUGLEVEL>=2)
-    { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); timer2(); }
+    { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); (void)timer2(); }
 
   /* Calcul de imax */
 
@@ -1563,12 +2218,12 @@ initzeta(GEN pol, long prec)
   zone  = switch_stack(NULL,i + 2*(N0+1) + 6*prec);
   zone1 = 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);
   tabcstni = (GEN*) cgetg(N0+1,t_VEC);
   p1 = ginv(cst);
   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 **********/
 
@@ -1595,7 +2250,7 @@ initzeta(GEN pol, long prec)
     if (k&1) ck_odd[k]=lneg((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[1] = serie_odd[1] = evalsigne(1)+evalvalp(1);
   i=0;
@@ -1656,7 +2311,8 @@ initzeta(GEN pol, long prec)
       t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1);
       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));
       }
       /* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */
@@ -1705,11 +2361,11 @@ initzeta(GEN pol, long prec)
     }
     /* use a parallel stack */
     z = i&1? zone1: zone0;
-    switch_stack(z, 1);
+    (void)switch_stack(z, 1);
     for (n=1; n<=N0; n++)
       if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]);
     /* come back */
-    switch_stack(z, 0);
+    (void)switch_stack(z, 0);
   }
   nfz[4] = (long) C;
   if (DEBUGLEVEL>=2) msgtimer("Cik");
@@ -1724,7 +2380,8 @@ gzetakall(GEN nfz, GEN s, long flag, long prec2)
   GEN resi,C,cst,cstlog,coeflog,cs,coef;
   GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2;
   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)
     err(talker,"not a zeta number field in zetakall");
@@ -1824,7 +2481,7 @@ gzetakall(GEN nfz, GEN s, long flag, long prec2)
 	for (k=1; k<=ru; 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);
 	}
 	val  = addis(val,1);