===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/buch4.c,v
retrieving revision 1.1.1.1
retrieving revision 1.2
diff -u -p -r1.1.1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/buch4.c	2001/10/02 11:17:03	1.1.1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/buch4.c	2002/09/11 07:26:50	1.2
@@ -1,4 +1,4 @@
-/* $Id: buch4.c,v 1.1.1.1 2001/10/02 11:17:03 noro Exp $
+/* $Id: buch4.c,v 1.2 2002/09/11 07:26:50 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -22,6 +22,11 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 #include "pari.h"
 #include "parinf.h"
 
+extern GEN to_polmod(GEN x, GEN mod);
+extern GEN hnfall0(GEN A, long remove);
+extern GEN get_theta_abstorel(GEN T, GEN pol, GEN k);
+extern GEN _rnfequation(GEN A, GEN B, long *pk, GEN *pLPRS);
+
 static long
 psquare(GEN a,GEN p)
 {
@@ -44,7 +49,8 @@ psquare(GEN a,GEN p)
 static long
 lemma6(GEN pol,GEN p,long nu,GEN x)
 {
-  long i,lambda,mu,ltop=avma;
+  long i, lambda, mu;
+  gpmem_t ltop=avma;
   GEN gx,gpx;
 
   for (i=lgef(pol)-2,gx=(GEN) pol[i+1]; i>1; i--)
@@ -65,7 +71,9 @@ lemma6(GEN pol,GEN p,long nu,GEN x)
 
 static long
 lemma7(GEN pol,long nu,GEN x)
-{ long i,odd4,lambda,mu,mnl,ltop=avma;
+{
+  long i,odd4,lambda,mu,mnl;
+  gpmem_t ltop=avma;
   GEN gx,gpx,oddgx;
 
   for (i=lgef(pol)-2,gx=(GEN) pol[i+1]; i>1; i--)
@@ -96,7 +104,8 @@ lemma7(GEN pol,long nu,GEN x)
 static long
 zpsol(GEN pol,GEN p,long nu,GEN pnu,GEN x0)
 {
-  long i,result,ltop=avma;
+  long i, result;
+  gpmem_t ltop=avma;
   GEN x,pnup;
 
   result = (cmpis(p,2)) ? lemma6(pol,p,nu,x0) : lemma7(pol,nu,x0);
@@ -137,7 +146,7 @@ qpsoluble(GEN pol,GEN p)
 static long
 psquarenf(GEN nf,GEN a,GEN pr)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long v;
   GEN norm;
 
@@ -169,7 +178,8 @@ check2(GEN nf, GEN a, GEN zinit)
 static long
 psquare2nf(GEN nf,GEN a,GEN pr,GEN zinit)
 {
-  long v, ltop = avma;
+  long v;
+  gpmem_t ltop = avma;
 
   if (gcmp0(a)) return 1;
   v = idealval(nf,a,pr); if (v&1) return 0;
@@ -182,7 +192,8 @@ psquare2nf(GEN nf,GEN a,GEN pr,GEN zinit)
 static long
 psquare2qnf(GEN nf,GEN a,GEN p,long q)
 {
-  long v, ltop=avma;
+  long v;
+  gpmem_t ltop=avma;
   GEN zinit = zidealstarinit(nf,idealpows(nf,p,q));
 
   v = check2(nf,a,zinit); avma = ltop; return v;
@@ -191,7 +202,8 @@ psquare2qnf(GEN nf,GEN a,GEN p,long q)
 static long
 lemma6nf(GEN nf,GEN pol,GEN p,long nu,GEN x)
 {
-  long i,lambda,mu,ltop=avma;
+  long i, lambda, mu;
+  gpmem_t ltop=avma;
   GEN gx,gpx;
 
   for (i=lgef(pol)-2,gx=(GEN) pol[i+1]; i>1; i--)
@@ -212,7 +224,8 @@ lemma6nf(GEN nf,GEN pol,GEN p,long nu,GEN x)
 static long
 lemma7nf(GEN nf,GEN pol,GEN p,long nu,GEN x,GEN zinit)
 {
-  long res,i,lambda,mu,q,ltop=avma;
+  long res, i, lambda, mu, q;
+  gpmem_t ltop=avma;
   GEN gx,gpx,p1;
 
   for (i=lgef(pol)-2, gx=(GEN) pol[i+1]; i>1; i--)
@@ -237,7 +250,7 @@ lemma7nf(GEN nf,GEN pol,GEN p,long nu,GEN x,GEN zinit)
     q=(nu<<1)-lambda; res=0;
   }
   if (q > itos((GEN) p[3])<<1) { avma=ltop; return -1; }
-  p1 = gmodulcp(gpuigs(gmul((GEN)nf[7],(GEN)p[2]), lambda), (GEN)nf[1]);
+  p1 = gmodulcp(gpowgs(gmul((GEN)nf[7],(GEN)p[2]), lambda), (GEN)nf[1]);
   if (!psquare2qnf(nf,gdiv(gx,p1), p,q)) res = -1;
   avma=ltop; return res;
 }
@@ -245,7 +258,8 @@ lemma7nf(GEN nf,GEN pol,GEN p,long nu,GEN x,GEN zinit)
 static long
 zpsolnf(GEN nf,GEN pol,GEN p,long nu,GEN pnu,GEN x0,GEN repr,GEN zinit)
 {
-  long i,result,ltop=avma;
+  long i, result;
+  gpmem_t ltop=avma;
   GEN pnup;
 
   nf=checknf(nf);
@@ -297,7 +311,7 @@ long
 qpsolublenf(GEN nf,GEN pol,GEN pr)
 {
   GEN repr,zinit,p1;
-  long ltop=avma;
+  gpmem_t ltop=avma;
 
   if (gcmp0(pol)) return 1;
   if (typ(pol)!=t_POL) err(notpoler,"qpsolublenf");
@@ -333,7 +347,7 @@ long
 zpsolublenf(GEN nf,GEN pol,GEN p)
 {
   GEN repr,zinit;
-  long ltop=avma;
+  gpmem_t ltop=avma;
 
   if (gcmp0(pol)) return 1;
   if (typ(pol)!=t_POL) err(notpoler,"zpsolublenf");
@@ -359,9 +373,13 @@ zpsolublenf(GEN nf,GEN pol,GEN p)
 static long
 hilb2nf(GEN nf,GEN a,GEN b,GEN p)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long rep;
-  GEN pol = coefs_to_pol(3, lift(a), zero, lift(b));
+  GEN pol;
+
+  if (typ(a) != t_POLMOD) a = basistoalg(nf, a);
+  if (typ(b) != t_POLMOD) b = basistoalg(nf, b);
+  pol = coefs_to_pol(3, lift(a), zero, lift(b));
   /* varn(nf.pol) = 0, pol is not a valid GEN  [as in Pol([x,x], x)].
    * But it is only used as a placeholder, hence it is not a problem */
 
@@ -373,9 +391,9 @@ hilb2nf(GEN nf,GEN a,GEN b,GEN p)
 long
 nfhilbertp(GEN nf,GEN a,GEN b,GEN pr)
 {
-  GEN ord, ordp, p, prhall,t;
+  GEN ord, ordp, T,p, modpr,t;
   long va, vb, rep;
-  ulong av = avma;
+  gpmem_t av = avma;
 
   if (gcmp0(a) || gcmp0(b)) err (talker,"0 argument in nfhilbertp");
   checkprimeid(pr); nf = checknf(nf);
@@ -394,9 +412,14 @@ nfhilbertp(GEN nf,GEN a,GEN b,GEN pr)
   /* quad. symbol is image of t by the quadratic character  */
   ord = subis( idealnorm(nf,pr), 1 ); /* |(O_K / pr)^*| */
   ordp= subis( p, 1);                 /* |F_p^*|        */
-  prhall = nfmodprinit(nf, pr);
-  t = element_powmodpr(nf, t, divii(ord, ordp), prhall); /* in F_p^* */
-  t = lift_intern((GEN)t[1]);
+  modpr = nf_to_ff_init(nf, &pr,&T,&p);
+  t = nf_to_ff(nf,t,modpr);
+  t = FpXQ_pow(t, diviiexact(ord, ordp), T,p); /* in F_p^* */
+  if (typ(t) == t_POL)
+  {
+    if (degpol(t)) err(bugparier,"nfhilbertp");
+    t = constant_term(t);
+  }
   rep = kronecker(t, p);
   avma = av; return rep;
 }
@@ -409,7 +432,7 @@ nfhilbertp(GEN nf,GEN a,GEN b,GEN pr)
 long
 nfhilbert(GEN nf,GEN a,GEN b)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long r1, i;
   GEN S, al, bl, ro;
 
@@ -457,7 +480,8 @@ nfhilbert0(GEN nf,GEN a,GEN b,GEN p)
 extern GEN isprincipalfact(GEN bnf,GEN P, GEN e, GEN C, long flag);
 extern GEN vconcat(GEN Q1, GEN Q2);
 extern GEN mathnfspec(GEN x, GEN *ptperm, GEN *ptdep, GEN *ptB, GEN *ptC);
-extern GEN factorback_i(GEN fa, GEN nf, int red);
+extern GEN factorback_i(GEN fa, GEN e, GEN nf, int red);
+extern GEN detcyc(GEN cyc);
 /* S a list of prime ideal in primedec format. Return res:
  * res[1] = generators of (S-units / units), as polynomials
  * res[2] = [perm, HB, den], for bnfissunit
@@ -469,10 +493,10 @@ extern GEN factorback_i(GEN fa, GEN nf, int red);
 GEN
 bnfsunit(GEN bnf,GEN S,long prec)
 {
-  ulong ltop = avma;
+  gpmem_t ltop = avma;
   long i,j,ls;
   GEN p1,nf,classgp,gen,M,U,H;
-  GEN sunit,card,sreg,res,pow,fa = cgetg(3, t_MAT);
+  GEN sunit,card,sreg,res,pow;
 
   if (typ(S) != t_VEC) err(typeer,"bnfsunit");
   bnf = checkbnf(bnf); nf=(GEN)bnf[7];
@@ -504,23 +528,19 @@ bnfsunit(GEN bnf,GEN S,long prec)
   card = gun;
   if (lg(H) > 1)
   { /* non trivial (rare!) */
-    GEN SNF, ClS = cgetg(4,t_VEC);
+    GEN D,U, ClS = cgetg(4,t_VEC);
 
-    SNF = smith2(H); p1 = (GEN)SNF[3];
-    card = dethnf_i(p1);
+    D = smithall(H, &U, NULL);
+    for(i=1; i<lg(D); i++)
+      if (gcmp1((GEN)D[i])) break;
+    setlg(D,i); D = mattodiagonal_i(D); /* cf smithrel */
+    card = detcyc(D);
     ClS[1] = (long)card; /* h */
-    for(i=1; i<lg(p1); i++)
-      if (gcmp1((GEN)p1[i])) break;
-    setlg(p1,i);
-    ClS[2]=(long)p1; /* cyc */
+    ClS[2] = (long)D; /* cyc */
 
-    p1=cgetg(i,t_VEC); pow=ZM_inv((GEN)SNF[1],gun);
-    fa[1] = (long)gen;
+    p1=cgetg(i,t_VEC); pow=ZM_inv(U,gun);
     for(i--; i; i--)
-    {
-      fa[2] = pow[i];
-      p1[i] = (long)factorback_i(fa, nf, 1);
-    }
+      p1[i] = (long)factorback_i(gen, (GEN)pow[i], nf, 1);
     ClS[3]=(long)p1; /* gen */
     res[5]=(long) ClS;
   }
@@ -546,7 +566,7 @@ bnfsunit(GEN bnf,GEN S,long prec)
     Sperm = cgetg(ls, t_VEC); sunit = cgetg(ls, t_VEC);
     for (i=1; i<ls; i++) Sperm[i] = S[perm[i]]; /* S o perm */
 
-    setlg(Sperm, lH); fa[1] = (long)Sperm;
+    setlg(Sperm, lH);
     for (i=1; i<lH; i++)
       sunit[i] = isprincipalfact(bnf,Sperm,(GEN)H[i],NULL,fl)[2];
     for (j=1; j<lB; j++,i++)
@@ -574,41 +594,26 @@ bnfsunit(GEN bnf,GEN S,long prec)
   return gerepilecopy(ltop,res);
 }
 
-/* cette fonction est l'equivalent de isunit, sauf qu'elle donne le resultat
- * avec des s-unites: si x n'est pas une s-unite alors issunit=[]~;
- * si x est une s-unite alors
- * x=\prod_{i=0}^r {e_i^issunit[i]}*prod{i=r+1}^{r+s} {s_i^issunit[i]}
- * ou les e_i sont les unites du corps (comme dans isunit)
- * et les s_i sont les s-unites calculees par sunit (dans le meme ordre).
- */
-GEN
-bnfissunit(GEN bnf,GEN suni,GEN x)
+static GEN
+make_unit(GEN bnf, GEN suni, GEN *px)
 {
-  long lB,cH,i,k,ls,tetpil, av = avma;
+  long lB, cH, i, k, ls;
   GEN den,gen,S,v,p1,xp,xm,xb,N,HB,perm;
 
-  bnf = checkbnf(bnf);
-  if (typ(suni)!=t_VEC || lg(suni)!=7) err(typeer,"bnfissunit");
-  switch (typ(x))
-  {
-    case t_INT: case t_FRAC: case t_FRACN:
-    case t_POL: case t_COL:
-      x = basistoalg(bnf,x); break;
-    case t_POLMOD: break;
-    default: err(typeer,"bnfissunit");
-  }
-  if (gcmp0(x)) return cgetg(1,t_COL);
-
+  if (gcmp0(*px)) return NULL;
   S = (GEN) suni[6]; ls=lg(S);
-  if (ls==1) return isunit(bnf,x);
+  if (ls==1) return cgetg(1, t_COL);
 
+  xb = algtobasis(bnf,*px); p1 = Q_denom(xb);
+  N = mulii(gnorm(gmul(*px,p1)), p1); /* relevant primes divide N */
+  if (is_pm1(N)) return zerocol(ls -1);
+
   p1 = (GEN)suni[2];
   perm = (GEN)p1[1];
-  HB = (GEN)p1[2]; den = (GEN)p1[3];
+  HB = (GEN)p1[2];
+  den = (GEN)p1[3];
   cH = lg(HB[1]) - 1;
   lB = lg(HB) - cH;
-  xb = algtobasis(bnf,x); p1 = denom(content(xb));
-  N = mulii(gnorm(gmul(x,p1)), p1); /* relevant primes divide N */
   v = cgetg(ls, t_VECSMALL);
   for (i=1; i<ls; i++)
   {
@@ -622,7 +627,7 @@ bnfissunit(GEN bnf,GEN suni,GEN x)
   for (i=1; i<=cH; i++)
   {
     GEN w = gdiv((GEN)v[i], den);
-    if (typ(w) != t_INT) { avma = av; return cgetg(1,t_COL); }
+    if (typ(w) != t_INT) return NULL;
     v[i] = (long)w;
   }
   p1 += cH;
@@ -634,174 +639,244 @@ bnfissunit(GEN bnf,GEN suni,GEN x)
   {
     k = -itos((GEN)v[i]); if (!k) continue;
     p1 = basistoalg(bnf, (GEN)gen[i]);
-    if (k > 0) xp = gmul(xp, gpuigs(p1, k));
-    else       xm = gmul(xm, gpuigs(p1,-k));
+    if (k > 0) xp = gmul(xp, gpowgs(p1, k));
+    else       xm = gmul(xm, gpowgs(p1,-k));
   }
-  if (xp != gun) x = gmul(x,xp);
-  if (xm != gun) x = gdiv(x,xm);
-  p1 = isunit(bnf,x);
-  if (lg(p1)==1) { avma = av; return cgetg(1,t_COL); }
-  tetpil=avma; return gerepile(av,tetpil,concat(p1,v));
+  if (xp != gun) *px = gmul(*px,xp);
+  if (xm != gun) *px = gdiv(*px,xm);
+  return v;
 }
 
-static void
-vecconcat(GEN bnf,GEN relnf,GEN vec,GEN *prod,GEN *S1,GEN *S2)
+/* cette fonction est l'equivalent de isunit, sauf qu'elle donne le resultat
+ * avec des s-unites: si x n'est pas une s-unite alors issunit=[]~;
+ * si x est une s-unite alors
+ * x=\prod_{i=0}^r {e_i^issunit[i]}*prod{i=r+1}^{r+s} {s_i^issunit[i]}
+ * ou les e_i sont les unites du corps (comme dans isunit)
+ * et les s_i sont les s-unites calculees par sunit (dans le meme ordre).
+ */
+GEN
+bnfissunit(GEN bnf,GEN suni,GEN x)
 {
-  long i;
+  gpmem_t av = avma;
+  GEN v, w;
 
-  for (i=1; i<lg(vec); i++)
-    if (signe(resii(*prod,(GEN)vec[i])))
-    {
-       *prod=mulii(*prod,(GEN)vec[i]);
-       *S1=concatsp(*S1,primedec(bnf,(GEN)vec[i]));
-       *S2=concatsp(*S2,primedec(relnf,(GEN)vec[i]));
-    }
+  bnf = checkbnf(bnf);
+  if (typ(suni)!=t_VEC || lg(suni)!=7) err(typeer,"bnfissunit");
+  switch (typ(x))
+  {
+    case t_INT: case t_FRAC: case t_FRACN:
+    case t_POL: case t_COL:
+      x = basistoalg(bnf,x); break;
+    case t_POLMOD: break;
+    default: err(typeer,"bnfissunit");
+  }
+  v = NULL;
+  if ( (w = make_unit(bnf, suni, &x)) ) v = isunit(bnf, x);
+  if (!v || lg(v) == 1) { avma = av; return cgetg(1,t_COL); }
+  return gerepileupto(av, concat(v, w));
 }
 
-/* bnf est le corps de base (buchinitfu).
- * ext definit l'extension relative:
- * ext[1] est une equation relative du corps,
- * telle qu'une de ses racines engendre le corps sur Q.
- * ext[2] exprime le generateur (y) du corps de base,
- * en fonction de la racine (x) de ext[1],
- * ext[3] est le buchinitfu (sur Q) de l'extension.
+static void
+pr_append(GEN nf, GEN rel, GEN p, GEN *prod, GEN *S1, GEN *S2)
+{
+  if (divise(*prod, p)) return;
+  *prod = mulii(*prod, p);
+  *S1 = concatsp(*S1, primedec(nf,p));
+  *S2 = concatsp(*S2, primedec(rel,p));
+}
 
- * si flag=0 c'est qu'on sait a l'avance que l'extension est galoisienne,
- * et dans ce cas la reponse est exacte.
- * si flag>0 alors on ajoue dans S tous les ideaux qui divisent p<=flag.
- * si flag<0 alors on ajoute dans S tous les ideaux qui divisent -flag.
+static void
+fa_pr_append(GEN nf,GEN rel,GEN N,GEN *prod,GEN *S1,GEN *S2)
+{
+  if (!is_pm1(N))
+  {
+    GEN v = (GEN)factor(N)[1];
+    long i, l = lg(v);
+    for (i=1; i<l; i++) pr_append(nf,rel,(GEN)v[i],prod,S1,S2);
+  }
+}
 
- * la reponse est un vecteur v a 2 composantes telles que
- * x=N(v[1])*v[2].
- * x est une norme ssi v[2]=1.
- */
+static GEN
+pol_up(GEN rnfeq, GEN x)
+{
+  long i, l = lgef(x);
+  GEN y = cgetg(l, t_POL); y[1] = x[1];
+  for (i=1; i<l; i++) y[i] = (long)rnfelementreltoabs(rnfeq, (GEN)x[i]);
+  return y;
+}
+
 GEN
-rnfisnorm(GEN bnf,GEN ext,GEN x,long flag,long PREC)
+rnfisnorminit(GEN T, GEN relpol, int galois)
 {
-  long lgsunitrelnf,i;
-  ulong ltop = avma;
-  GEN relnf,aux,vec,tors,xnf,H,Y,M,A,suni,sunitrelnf,sunitnormnf,prod;
-  GEN res = cgetg(3,t_VEC), S1,S2;
+  gpmem_t av = avma;
+  long i, l, drel;
+  GEN prod, S1, S2, gen, cyc, bnf, nf, nfrel, rnfeq, rel, res, k, polabs;
+  GEN y = cgetg(9, t_VEC);
 
-  if (typ(ext)!=t_VEC || lg(ext)!=4) err (typeer,"bnfisnorm");
-  if (typ(x)!=t_POL) x = basistoalg(bnf,x);
-  bnf = checkbnf(bnf); relnf = (GEN)ext[3];
-  if (gcmp0(x) || gcmp1(x) || (gcmp_1(x) && (degpol(ext[1])&1)))
-  {
-    avma = (long)res; res[1]=lcopy(x); res[2]=un; return res;
-  }
+  T = get_bnfpol(T, &bnf, &nf);
+  if (!bnf) bnf = bnfinit0(nf? nf: T, 1, NULL, DEFAULTPREC);
+  if (!nf) nf = checknf(bnf);
 
-/* construction de l'ensemble S des ideaux
-   qui interviennent dans les solutions */
+  relpol = get_bnfpol(relpol, &rel, &nfrel);
+  drel = degpol(relpol);
 
-  prod=gun; S1=S2=cgetg(1,t_VEC);
-  if (!gcmp1(gmael3(relnf,8,1,1)))
+  rnfeq = NULL; /* no reltoabs needed */
+  if (degpol(nf[1]) == 1)
+  { /* over Q */
+    polabs = lift(relpol);
+    k = gzero;
+  }
+  else
   {
-    GEN genclass=gmael3(relnf,8,1,3);
-    vec=cgetg(1,t_VEC);
-    for(i=1;i<lg(genclass);i++)
-      if (!gcmp1(ggcd(gmael4(relnf,8,1,2,i), stoi(degpol(ext[1])))))
-        vec=concatsp(vec,(GEN)factor(gmael3(genclass,i,1,1))[1]);
-    vecconcat(bnf,relnf,vec,&prod,&S1,&S2);
+    if (galois == 2 && drel > 2)
+    { /* needs reltoabs */
+      rnfeq = rnfequation2(bnf, relpol);
+      polabs = (GEN)rnfeq[1];
+      k =      (GEN)rnfeq[3];
+    }
+    else
+    {
+      long sk;
+      polabs = _rnfequation(bnf, relpol, &sk, NULL);
+      k = stoi(sk);
+    }
   }
+  if (!rel || !gcmp0(k)) rel = bnfinit0(polabs, 1, NULL, nfgetprec(nf));
+  if (!nfrel) nfrel = checknf(rel);
 
-  if (flag>1)
+  if (galois < 0 || galois > 2) err(flagerr, "rnfisnorminit");
+  if (galois == 2)
   {
-    for (i=2; i<=flag; i++)
-      if (isprime(stoi(i)) && signe(resis(prod,i)))
-      {
-	prod=mulis(prod,i);
-	S1=concatsp(S1,primedec(bnf,stoi(i)));
-	S2=concatsp(S2,primedec(relnf,stoi(i)));
-      }
+    GEN P = rnfeq? pol_up(rnfeq, relpol): relpol;
+    galois = nfisgalois(gsubst(nfrel, varn(P), polx[varn(T)]), P);
   }
-  else if (flag<0)
-    vecconcat(bnf,relnf,(GEN)factor(stoi(-flag))[1],&prod,&S1,&S2);
 
-  if (flag)
+  prod = gun; S1 = S2 = cgetg(1, t_VEC);
+  res = gmael(rel,8,1);
+  cyc = (GEN)res[2];
+  gen = (GEN)res[3]; l = lg(cyc);
+  for(i=1; i<l; i++)
   {
-    GEN normdiscrel=divii(gmael(relnf,7,3),
-			  gpuigs(gmael(bnf,7,3),lg(ext[1])-3));
-    vecconcat(bnf,relnf,(GEN) factor(absi(normdiscrel))[1],
-	      &prod,&S1,&S2);
+    if (cgcd(smodis((GEN)cyc[i], drel), drel) == 1) break;
+    fa_pr_append(nf,rel,gmael3(gen,i,1,1),&prod,&S1,&S2);
   }
-  vec=(GEN) idealfactor(bnf,x)[1]; aux=cgetg(2,t_VEC);
-  for (i=1; i<lg(vec); i++)
-    if (signe(resii(prod,gmael(vec,i,1))))
-    {
-      aux[1]=vec[i]; S1=concatsp(S1,aux);
-    }
-  xnf=lift(x);
-  xnf=gsubst(xnf,varn(xnf),(GEN)ext[2]);
-  vec=(GEN) idealfactor(relnf,xnf)[1];
-  for (i=1; i<lg(vec); i++)
-    if (signe(resii(prod,gmael(vec,i,1))))
-    {
-      aux[1]=vec[i]; S2=concatsp(S2,aux);
-    }
+  if (!galois)
+  {
+    GEN Ndiscrel = diviiexact((GEN)nfrel[3], gpowgs((GEN)nf[3], drel));
+    fa_pr_append(nf,rel,absi(Ndiscrel), &prod,&S1,&S2);
+  }
 
-  res[1]=un; res[2]=(long)x;
-  tors=cgetg(2,t_VEC); tors[1]=mael3(relnf,8,4,2);
+  y[1] = (long)bnf;
+  y[2] = (long)rel;
+  y[3] = (long)relpol;
+  y[4] = (long)get_theta_abstorel(T, relpol, k);
+  y[5] = (long)prod;
+  y[6] = (long)S1;
+  y[7] = (long)S2;
+  y[8] = lstoi(galois); return gerepilecopy(av, y);
+}
 
-  /* calcul sur les S-unites */
+/* T as output by rnfisnorminit
+ * if flag=0 assume extension is Galois (==> answer is unconditional)
+ * if flag>0 add to S all primes dividing p <= flag
+ * if flag<0 add to S all primes dividing abs(flag)
 
-  suni=bnfsunit(bnf,S1,PREC);
-  A=lift(bnfissunit(bnf,suni,x));
-  sunitrelnf=(GEN) bnfsunit(relnf,S2,PREC)[1];
-  if (lg(sunitrelnf)>1)
+ * answer is a vector v = [a,b] such that
+ * x = N(a)*b and x is a norm iff b = 1  [assuming S large enough] */
+GEN
+rnfisnorm(GEN T, GEN x, long flag)
+{
+  gpmem_t av = avma;
+  GEN bnf = (GEN)T[1], rel = (GEN)T[2], relpol = (GEN)T[3], theta = (GEN)T[4];
+  GEN nf, aux, H, Y, M, A, suni, sunitrel, futu, tu, w;
+  GEN prod, S1, S2;
+  GEN res = cgetg(3,t_VEC);
+  long L, i, drel, itu;
+
+  if (typ(T) != t_VEC || lg(T) != 9)
+    err(talker,"please apply rnfisnorminit first");
+  bnf = checkbnf(bnf);
+  rel = checkbnf(rel);
+  nf = checknf(bnf);
+  x = basistoalg(nf,x);
+  if (typ(x) != t_POLMOD) err(typeer, "rnfisnorm");
+  drel = degpol(relpol);
+  if (gcmp0(x) || gcmp1(x) || (gcmp_1(x) && odd(drel)))
   {
-    sunitrelnf=lift(basistoalg(relnf,sunitrelnf));
-    sunitrelnf=concatsp(tors,sunitrelnf);
+    res[1] = (long)x;
+    res[2] = un; return res;
   }
-  else sunitrelnf=tors;
-  aux=(GEN)relnf[8];
-  if (lg(aux)>=6) aux=(GEN)aux[5];
-  else
+
+  /* build set T of ideals involved in the solutions */
+  prod = (GEN)T[5];
+  S1   = (GEN)T[6];
+  S2   = (GEN)T[7];
+  if (flag && !gcmp0((GEN)T[8]))
+    err(warner,"useless flag in rnfisnorm: the extension is Galois");
+  if (flag > 0)
   {
-    aux=buchfu(relnf);
-    if(gcmp0((GEN)aux[2]))
-      err(precer,"bnfisnorm, please increase precision and try again");
-    aux=(GEN)aux[1];
+    byteptr d = diffptr;
+    long p = 0;
+    if (maxprime() < flag) err(primer1);
+    for(;;)
+    {
+      NEXT_PRIME_VIADIFF(p, d);
+      if (p > flag) break;
+      pr_append(nf,rel,stoi(p),&prod,&S1,&S2);
+    }
   }
-  if (lg(aux)>1)
-    sunitrelnf=concatsp(aux,sunitrelnf);
-  lgsunitrelnf=lg(sunitrelnf);
-  M=cgetg(lgsunitrelnf+1,t_MAT);
-  sunitnormnf=cgetg(lgsunitrelnf,t_VEC);
-  for (i=1; i<lgsunitrelnf; i++)
+  else if (flag < 0)
+    fa_pr_append(nf,rel,stoi(-flag),&prod,&S1,&S2);
+  /* overkill: prime ideals dividing x would be enough */
+  fa_pr_append(nf,rel,idealnorm(nf,x), &prod,&S1,&S2);
+
+  /* computation on T-units */
+  w  = gmael3(rel,8,4,1);
+  tu = gmael3(rel,8,4,2);
+  futu = concatsp(check_units(rel,"rnfisnorm"), tu);
+  suni = bnfsunit(bnf,S1,3);
+  sunitrel = (GEN)bnfsunit(rel,S2,3)[1];
+  if (lg(sunitrel) > 1) sunitrel = lift_intern(basistoalg(rel,sunitrel));
+  sunitrel = concatsp(futu, sunitrel);
+
+  A = lift(bnfissunit(bnf,suni,x));
+  L = lg(sunitrel);
+  itu = lg(nf[6])-1; /* index of torsion unit in bnfsunit(nf) output */
+  M = cgetg(L+1,t_MAT);
+  for (i=1; i<L; i++)
   {
-    sunitnormnf[i]=lnorm(gmodulcp((GEN) sunitrelnf[i],(GEN)ext[1]));
-    M[i]=llift(bnfissunit(bnf,suni,(GEN) sunitnormnf[i]));
+    GEN u = poleval((GEN)sunitrel[i], theta); /* abstorel */
+    if (typ(u) != t_POLMOD) u = to_polmod(u, (GEN)theta[1]);
+    sunitrel[i] = (long)u;
+    u = bnfissunit(bnf,suni, gnorm(u));
+    if (lg(u) == 1) err(bugparier,"rnfisnorm");
+    u[itu] = llift((GEN)u[itu]); /* lift root of 1 part */
+    M[i] = (long)u;
   }
-  M[lgsunitrelnf]=lgetg(lg(A),t_COL);
-  for (i=1; i<lg(A); i++) mael(M,lgsunitrelnf,i)=zero;
-  mael(M,lgsunitrelnf,lg(mael(bnf,7,6))-1)=mael3(bnf,8,4,1);
-  H=hnfall(M); Y=inverseimage(gmul(M,(GEN) H[2]),A);
-  Y=gmul((GEN) H[2],Y);
-  for (aux=(GEN)res[1],i=1; i<lgsunitrelnf; i++)
-    aux=gmul(aux,gpuigs(gmodulcp((GEN) sunitrelnf[i],(GEN)ext[1]),
-                        itos(gfloor((GEN)Y[i]))));
-  x = gdiv(x,gnorm(gmodulcp(lift(aux),(GEN)ext[1])));
-  if (typ(x) == t_POLMOD && (typ(x[2]) != t_POL || lgef(x[2]) == 3))
+  aux = zerocol(lg(A)-1); aux[itu] = (long)w;
+  M[L] = (long)aux;
+  H = hnfall0(M, 0);
+  Y = gmul((GEN)H[2], inverseimage((GEN)H[1],A));
+  /* Y: sols of MY = A over Q */
+  setlg(Y, L);
+  aux = factorback(sunitrel, gfloor(Y));
+  x = gdiv(x, gnorm(gmodulcp(lift_intern(aux), relpol)));
+  if (typ(x) == t_POLMOD && (typ(x[2]) != t_POL || !degpol(x[2])))
   {
     x = (GEN)x[2]; /* rational number */
     if (typ(x) == t_POL) x = (GEN)x[2];
   }
-  res[1]=(long)aux;
-  res[2]=(long)x;
-  return gerepilecopy(ltop,res);
+  if (typ(aux) == t_POLMOD && degpol(nf[1]) == 1)
+    aux[2] = (long)lift_intern((GEN)aux[2]);
+  res[1] = (long)aux;
+  res[2] = (long)x;
+  return gerepilecopy(av, res);
 }
 
 GEN
 bnfisnorm(GEN bnf,GEN x,long flag,long PREC)
 {
-  long ltop = avma, lbot;
-  GEN ext = cgetg(4,t_VEC);
-
-  bnf = checkbnf(bnf);
-  ext[1] = mael(bnf,7,1);
-  ext[2] = zero;
-  ext[3] = (long) bnf;
-  bnf = buchinitfu(polx[MAXVARN],NULL,NULL,0); lbot = avma;
-  return gerepile(ltop,lbot,rnfisnorm(bnf,ext,x,flag,PREC));
+  gpmem_t av = avma;
+  GEN T = rnfisnorminit(polx[MAXVARN], bnf, flag == 0? 1: 2);
+  return gerepileupto(av, rnfisnorm(T, x, flag));
 }