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

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

version 1.1, 2001/10/02 11:17:11 version 1.2, 2002/09/11 07:27:05
Line 19  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
Line 19  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /*                                                                 */  /*                                                                 */
 /*******************************************************************/  /*******************************************************************/
 #include "pari.h"  #include "pari.h"
   #include "parinf.h"
   extern GEN gmul_mati_smallvec(GEN x, GEN y);
 extern GEN check_and_build_cycgen(GEN bnf);  extern GEN check_and_build_cycgen(GEN bnf);
 extern GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);  extern GEN get_arch_real(GEN nf,GEN x,GEN *emb,long prec);
 extern GEN vecmul(GEN x, GEN y);  
 extern GEN vecinv(GEN x);  
 extern GEN T2_from_embed_norm(GEN x, long r1);  extern GEN T2_from_embed_norm(GEN x, long r1);
 extern GEN pol_to_vec(GEN x, long N);  extern GEN vconcat(GEN A, GEN B);
 extern GEN famat_to_nf(GEN nf, GEN f);  extern long int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx);
   
 static long rc,ell,degK,degKz,m,d,vnf,dv;  extern GEN famat_inv(GEN f);
 static GEN matexpoteta1,nf,raycyc,polnf;  extern GEN famat_pow(GEN f, GEN n);
 static GEN bnfz,nfz,U,uu,gell,cyc,gencyc,vecalpha,R,g;  extern GEN famat_mul(GEN f, GEN g);
 static GEN listmod,listprSp,listbid,listunif,listellrank;  extern GEN famat_reduce(GEN fa);
 static GEN listbidsup,listellranksup,vecw;  extern GEN famat_ideallog(GEN nf, GEN g, GEN e, GEN bid);
   extern GEN to_famat(GEN g, GEN e);
   
 /* row vector B x matrix T : c_j=prod_i (b_i^T_ij) */  typedef struct {
 static GEN    GEN x;  /* tau ( Mod(x, nf.pol) ) */
 groupproduct(GEN B, GEN T)    GEN zk; /* action of tau on nf.zk (as t_MAT) */
 {  } tau_s;
   long lB,lT,i,j;  
   GEN c,p1;  
   
   lB=lg(B)-1;  typedef struct {
   lT=lg(T)-1;    GEN polnf, invexpoteta1;
   c=cgetg(lT+1,t_VEC);    tau_s *tau;
   for (j=1; j<=lT; j++)    long m;
   {  } toK_s;
     p1=gun;  
     for (i=1; i<=lB; i++) p1=gmul(p1,gpui((GEN)B[i],gcoeff(T,i,j),0));  
     c[j]=(long)p1;  
   }  
   return c;  
 }  
   
 static GEN  long
 grouppows(GEN B, long ex)  prank(GEN cyc, long ell)
 {  {
   long lB = lg(B),j;    long i;
   GEN c;    for (i=1; i<lg(cyc); i++)
       if (smodis((GEN)cyc[i],ell)) break;
   c = cgetg(lB,t_VEC);    return i-1;
   for (j=1; j<lB; j++) c[j] = lpowgs((GEN)B[j], ex);  
   return c;  
 }  }
   
   /* increment y, which runs through [0,ell-1]^(k-1). Return 0 when done. */
 static int  static int
 ok_x(GEN X, GEN arch, GEN vecmunit2, GEN msign)  increment(GEN y, long k, long ell)
 {  {
   long i, l = lg(vecmunit2);    long i = k, j;
   GEN p1;    do
   for (i=1; i<l; i++)  
   {    {
     p1 = FpV_red(gmul((GEN)vecmunit2[i], X), gell);      if (--i == 0) return 0;
     if (gcmp0(p1)) return 0;      y[i]++;
   }    } while (y[i] >= ell);
   p1 = lift(gmul(msign,X)); settyp(p1,t_VEC);    for (j = i+1; j < k; j++) y[j] = 0;
   return gegal(p1, arch);    return 1;
 }  }
   
 static GEN  /* as above, y increasing (y[i] <= y[i+1]) */
 get_listx(GEN arch,GEN msign,GEN munit,GEN vecmunit2,long ell,long lSp,long nbvunit)  static int
   increment_inc(GEN y, long k, long ell)
 {  {
   GEN kermat,p2,X,y, listx = cgetg(1,t_VEC);    long i = k, j;
   long i,j,cmpt,lker;    do
   
   kermat = FpM_ker(munit,gell); lker=lg(kermat)-1;  
   if (!lker) return listx;  
   y = cgetg(lker,t_VECSMALL);  
   for (i=1; i<lker; i++) y[i] = 0;  
   for(;;)  
   {    {
     p2 = cgetg(2,t_VEC);      if (--i == 0) return 0;
     X = (GEN)kermat[lker];      y[i]++;
     for (j=1; j<lker; j++) X = gadd(X, gmulsg(y[j],(GEN)kermat[j]));    } while (y[i] >= ell);
     X = FpV_red(X, gell);    for (j = i+1; j < k; j++) y[j] = y[i];
     cmpt = 0; for (j=1; j<=lSp; j++) if (gcmp0((GEN)X[j+nbvunit])) cmpt++;    return 1;
     if (!cmpt)  
     {  
       if (ok_x(X, arch, vecmunit2, msign))  
         { p2[1] = (long)X; listx = concatsp(listx,p2); }  
     }  
     if (!lSp)  
     {  
       j = 1; while (j<lker && !y[j]) j++;  
       if (j<lker && y[j] == 1)  
       {  
         X = gsub(X,(GEN)kermat[lker]);  
         if (ok_x(X, arch, vecmunit2, msign))  
           { p2[1] = (long)X; listx = concatsp(listx,p2); }  
       }  
     }  
     i = lker;  
     do  
     {  
       i--; if (!i) return listx;  
       if (i < lker-1) y[i+1] = 0;  
       y[i]++;  
     }  
     while (y[i] >= ell);  
   }  
 }  }
   
 static GEN  static int
 reducealpha(GEN nf,GEN x,GEN gell)  ok_congruence(GEN X, GEN ell, long lW, GEN vecMsup)
 /* etant donne un nombre algebrique x du corps nf -- non necessairement  
    entier -- calcule un entier algebrique de la forme x*g^gell */  
 {  {
   long tx=typ(x),i;    long i, l;
   GEN den,fa,fac,ep,p1,y;    if (gcmp0(X)) return 0;
     l = lg(X);
   nf = checknf(nf);    for (i=lW; i<l; i++)
   if (tx==t_POL || tx==t_POLMOD) y = algtobasis(nf,x);      if (gcmp0((GEN)X[i])) return 0;
   else { y = x; x = basistoalg(nf,y); }    l = lg(vecMsup);
   den = denom(y);    for (i=1; i<l; i++)
   if (gcmp1(den)) return x;      if (gcmp0(FpV_red(gmul((GEN)vecMsup[i],X), ell))) return 0;
   fa = decomp(den); fac = (GEN)fa[1];ep = (GEN)fa[2];    return 1;
   p1 = gun;  
   for (i=1; i<lg(fac); i++)  
     p1 = mulii(p1, powgi((GEN)fac[i], gceil(gdiv((GEN)ep[i],gell))));  
   return gmul(powgi(p1, gell), x);  
 }  }
   
 long  static int
 ellrank(GEN cyc, long ell)  ok_sign(GEN X, GEN msign, GEN arch)
 {  {
   long i;    GEN p1 = lift(gmul(msign,X)); settyp(p1,t_VEC);
   for (i=1; i<lg(cyc); i++)    return gegal(p1, arch);
     if (smodis((GEN)cyc[i],ell)) break;  
   return i-1;  
 }  }
   
   /* REDUCTION MOD ell-TH POWERS */
   
 static GEN  static GEN
 no_sol(long all, int i)  fix_be(GEN bnfz, GEN be, GEN u)
 {  {
   if (!all) err(talker,"bug%d in kummer",i);    GEN nf = checknf(bnfz), fu = gmael(bnfz,8,5);
   return cgetg(1,t_VEC);    return element_mul(nf, be, factorbackelt(fu, u, nf));
 }  }
   
 static void  static GEN
 _append(GEN x, GEN t)  logarch2arch(GEN x, long r1, long prec)
 {  {
   long l = lg(x); x[l] = (long)t; setlg(x, l+1);    long i, lx = lg(x), tx = typ(x);
     GEN y = cgetg(lx, tx);
   
     if (tx == t_MAT)
     {
       for (i=1; i<lx; i++) y[i] = (long)logarch2arch((GEN)x[i], r1, prec);
       return y;
     }
     for (i=1; i<=r1;i++) y[i] = lexp((GEN)x[i],prec);
     for (   ; i<lx; i++) y[i] = lexp(gmul2n((GEN)x[i],-1),prec);
     return y;
 }  }
   
 /* si all!=0, donne toutes les equations correspondant a un sousgroupe  /* multiply be by ell-th powers of units as to find small L2-norm for new be */
    de determinant all (i.e. de degre all) */  
 static GEN  static GEN
 rnfkummersimple(GEN bnr, GEN subgroup, long all)  reducebetanaive(GEN bnfz, GEN be, GEN b, GEN ell)
 {  {
   long r1,degnf,ell,j,i,l;    long i,k,n,ru,r1, prec = nfgetprec(bnfz);
   long nbgenclK,lSml2,lSl,lSp,rc,nbvunit;    GEN z,p1,p2,nmax,c, nf = checknf(bnfz);
   long lastbid,nbcol,nblig,k,llistx,e,vp,vd;  
   
   GEN bnf,nf,classgroup,cyclic,bid,ideal,arch,cycgen,gell,p1,p2,p3;    r1 = nf_get_r1(nf);
   GEN cyclicK,genK,listgamma,listalpha,fa,prm,expom,Sm,Sml1,Sml2,Sl,ESml2;    if (!b)
   GEN p,factell,Sp,vecbeta,matexpo,vunit,id,pr,vsml2,vecalpha0;  
   GEN cycpro,munit,vecmunit2,msign,archartif,listx,listal,listg,listgamma0;  
   GEN vecbeta0;  
   
   checkbnrgen(bnr);  
   bnf = (GEN)bnr[1];  
   nf = (GEN)bnf[7]; r1 = nf_get_r1(nf);  
   polnf = (GEN)nf[1]; vnf = varn(polnf); degnf = degpol(polnf);  
   if (vnf==0) err(talker,"main variable in kummer must not be x");  
   
   p1 = conductor(bnr,all ? gzero : subgroup,2);  
   bnr = (GEN)p1[2];  
   if (!all) subgroup = (GEN)p1[3];  
   classgroup = (GEN)bnr[5];  
   cyclic = (GEN)classgroup[2];  
   bid = (GEN)bnr[2];  
   ideal= gmael(bid,1,1);  
   arch = gmael(bid,1,2); /* this is the conductor */  
   if (!all && gcmp0(subgroup)) subgroup = diagonal(cyclic);  
   gell = all? stoi(all): det(subgroup);  
   ell = itos(gell);  
   
   cyclicK= gmael3(bnf,8,1,2); rc = ellrank(cyclicK,ell);  
   genK   = gmael3(bnf,8,1,3); nbgenclK = lg(genK)-1;  
   listgamma0=cgetg(nbgenclK+1,t_VEC);  
   listgamma=cgetg(nbgenclK+1,t_VEC);  
   vecalpha0=cgetg(rc+1,t_VEC);  
   listalpha=cgetg(rc+1,t_VEC);  
   cycgen = check_and_build_cycgen(bnf);  
   p1 = gmul(gell,ideal);  
   for (i=1; i<=rc; i++)  
   {    {
     p3 = basistoalg(nf,idealcoprime(nf,(GEN)genK[i],p1));      if (typ(be) != t_COL) be = algtobasis(nf, be);
     p2 = basistoalg(nf, famat_to_nf(nf, (GEN)cycgen[i]));      b = gmul(gmael(nf,5,1), be);
     listgamma[i] = listgamma0[i] = linv(p3);  
     vecalpha0[i] = (long)p2;  
     listalpha[i] = lmul((GEN)vecalpha0[i], powgi(p3,(GEN)cyclicK[i]));  
   }    }
   for (   ; i<=nbgenclK; i++)    n = max((itos(ell)>>1), 3);
     z = cgetg(n+1, t_VEC);
     c = gmul(greal((GEN)bnfz[3]), ell);
     c = logarch2arch(c, r1, prec); /* = embeddings of fu^ell */
     c = gprec_w(gnorm(c), DEFAULTPREC);
     b = gprec_w(gnorm(b), DEFAULTPREC); /* need little precision */
     z[1] = (long)concatsp(c, vecinv(c));
     for (k=2; k<=n; k++) z[k] = (long) vecmul((GEN)z[1], (GEN)z[k-1]);
     nmax = T2_from_embed_norm(b, r1);
     ru = lg(c)-1; c = zerovec(ru);
     for(;;)
   {    {
     long k;      GEN B = NULL;
     p3 = basistoalg(nf,idealcoprime(nf,(GEN)genK[i],p1));      long besti = 0, bestk = 0;
     p2 = basistoalg(nf, famat_to_nf(nf, (GEN)cycgen[i]));      for (k=1; k<=n; k++)
     k = itos(mpinvmod((GEN)cyclicK[i], gell));        for (i=1; i<=ru; i++)
     p2 = gpowgs(p2,k);  
     listgamma0[i]= (long)p2;  
     listgamma[i] = lmul(p2, gpowgs(p3, k * itos((GEN)cyclicK[i]) - 1));  
   }  
   fa = (GEN)bid[3];  
   prm   = (GEN)fa[1];  
   expom = (GEN)fa[2]; l = lg(prm);  
   Sm  = cgetg(l,t_VEC); setlg(Sm,1);  
   Sml1= cgetg(l,t_VEC); setlg(Sml1,1);  
   Sml2= cgetg(l,t_VEC); setlg(Sml2,1);  
   Sl  = cgetg(l+degnf,t_VEC); setlg(Sl,1);  
   ESml2=cgetg(l,t_VECSMALL); setlg(ESml2,1);  
   for (i=1; i<l; i++)  
   {  
     pr = (GEN)prm[i]; p = (GEN)pr[1]; e = itos((GEN)pr[3]);  
     vp = itos((GEN)expom[i]);  
     if (!egalii(p,gell))  
     {  
       if (vp != 1) return no_sol(all,1);  
       _append(Sm, pr);  
     }  
     else  
     {  
       vd = (vp-1)*(ell-1)-ell*e;  
       if (vd > 0) return no_sol(all,4);  
       if (vd==0)  
         _append(Sml1, pr);  
       else  
       {        {
         if (vp==1) return no_sol(all,2);          p1 = vecmul(b, gmael(z,k,i));    p2 = T2_from_embed_norm(p1,r1);
         _append(Sml2, pr);          if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk = k; continue; }
         _append(ESml2,(GEN)vp);          p1 = vecmul(b, gmael(z,k,i+ru)); p2 = T2_from_embed_norm(p1,r1);
           if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk =-k; }
       }        }
     }      if (!B) break;
       b = B; c[besti] = laddis((GEN)c[besti], bestk);
   }    }
   factell = primedec(nf,gell); l = lg(factell);    if (DEBUGLEVEL) fprintferr("naive reduction mod U^l: unit exp. = %Z\n",c);
   for (i=1; i<l; i++)    return fix_be(bnfz, be, gmul(ell,c));
   }
   
   static GEN
   reduce_mod_Qell(GEN bnfz, GEN be, GEN gell)
   {
     GEN c, fa;
     if (typ(be) != t_COL) be = algtobasis(bnfz, be);
     be = primitive_part(be, &c);
     if (c)
   {    {
     pr = (GEN)factell[i];      fa = factor(c);
     if (!idealval(nf,ideal,pr)) _append(Sl, pr);      fa[2] = (long)FpV_red((GEN)fa[2], gell);
       c = factorback(fa, NULL);
       be = gmul(be, c);
   }    }
   lSml2=lg(Sml2)-1; lSl=lg(Sl)-1;    return be;
   Sp = concatsp(Sm,Sml1); lSp = lg(Sp)-1;  }
   vecbeta = cgetg(lSp+1,t_VEC); matexpo=cgetg(lSp+1,t_MAT);  
   vecbeta0= cgetg(lSp+1,t_VEC);  /* return q, q^n r = x, v_pr(r) < n for all pr. Insist q is a genuine n-th
   for (j=1; j<=lSp; j++)   * root (i.e r = 1) if strict != 0. */
   static GEN
   idealsqrtn(GEN nf, GEN x, GEN gn, int strict)
   {
     long i, l, n = itos(gn);
     GEN fa, q, Ex, Pr;
   
     fa = idealfactor(nf, x);
     Pr = (GEN)fa[1]; l = lg(Pr);
     Ex = (GEN)fa[2]; q = NULL;
     for (i=1; i<l; i++)
   {    {
     p1=isprincipalgenforce(bnf,(GEN)Sp[j]);      long ex = itos((GEN)Ex[i]);
     p2=basistoalg(nf,(GEN)p1[2]);      GEN e = stoi(ex / n);
     for (i=1; i<=rc; i++)      if (strict && ex % n) err(talker,"not an n-th power in idealsqrtn");
       p2 = gmul(p2, powgi((GEN)listgamma[i],gmael(p1,1,i)));      if (q) q = idealmulpowprime(nf, q, (GEN)Pr[i], e);
     p3 = p2;      else   q = idealpow(nf, (GEN)Pr[i], e);
     for (   ; i<=nbgenclK; i++)  
     {  
       p2 = gmul(p2, powgi((GEN)listgamma[i],gmael(p1,1,i)));  
       p3 = gmul(p3, powgi((GEN)listgamma0[i],gmael(p1,1,i)));  
     }  
     matexpo[j]=(long)p1[1];  
     vecbeta[j]=(long)p2; /* attention, ceci sont les beta modifies */  
     vecbeta0[j]=(long)p3;  
   }    }
   listg = gmodulcp(concatsp(gmael(bnf,8,5),gmael3(bnf,8,4,2)),polnf);    return q? q: gun;
   vunit = concatsp(listg, listalpha);  }
   nbvunit=lg(vunit)-1;  
   id=idmat(degnf);  static GEN
   for (i=1; i<=lSl; i++)  reducebeta(GEN bnfz, GEN be, GEN ell)
   {
     long j,ru, prec = nfgetprec(bnfz);
     GEN emb,z,u,matunit, nf = checknf(bnfz);
   
     if (DEBUGLEVEL>1) fprintferr("reducing beta = %Z\n",be);
     /* reduce mod Q^ell */
     be = reduce_mod_Qell(nf, be, ell);
     /* reduce l-th root */
     z = idealsqrtn(nf, be, ell, 0);
     if (typ(z) == t_MAT && !gcmp1(gcoeff(z,1,1)))
   {    {
     pr = (GEN)Sl[i]; e = itos((GEN)pr[3]);      z = idealred_elt(nf, z);
     id = idealmul(nf, id, idealpows(nf,pr,(ell*e)/(ell-1)));      be = element_div(nf, be, element_pow(nf, z, ell));
       /* make be integral */
       be = reduce_mod_Qell(nf, be, ell);
   }    }
   vsml2=cgetg(lSml2+1,t_VEC);    if (DEBUGLEVEL>1) fprintferr("beta reduced via ell-th root = %Z\n",be);
   for (i=1; i<=lSml2; i++)  
     matunit = gmul(greal((GEN)bnfz[3]), ell); /* log. embeddings of fu^ell */
     for (;;)
   {    {
     pr = (GEN)Sml2[i]; e = itos((GEN)pr[3]);      z = get_arch_real(nf, be, &emb, prec);
     p1=idealpows(nf,pr, (e*ell)/(ell-1) + 1 - ESml2[i]);      if (z) break;
     id = idealmul(nf,id,p1);      prec = (prec-1)<<1;
     p2 = idealmul(nf,p1,pr);      if (DEBUGLEVEL) err(warnprec,"reducebeta",prec);
     p3 = zidealstarinitgen(nf,p2);      nf = nfnewprec(nf,prec);
     vsml2[i] = (long)p3;  
     cycpro = gmael(p3,2,2); l = lg(cycpro)-1;  
     if (! gdivise((GEN)cycpro[l],gell)) err(talker,"bug5 in kummer");  
   }    }
   bid = zidealstarinitgen(nf,id);    z = concatsp(matunit, z);
   lastbid = ellrank(gmael(bid,2,2), ell);    u = lllintern(z, 100, 1, prec);
   nbcol=nbvunit+lSp; nblig=lastbid+rc;    if (u)
   munit=cgetg(nbcol+1,t_MAT); vecmunit2=cgetg(lSml2+1,t_VEC);  
   msign=cgetg(nbcol+1,t_MAT);  
   for (k=1; k<=lSml2; k++) vecmunit2[k]=lgetg(nbcol+1,t_MAT);  
   archartif=cgetg(r1+1,t_VEC); for (j=1; j<=r1; j++) archartif[j]=un;  
   for (j=1; j<=nbvunit; j++)  
   {    {
     p1 = zideallog(nf,(GEN)vunit[j],bid);      ru = lg(u);
     p2 = cgetg(nblig+1,t_COL); munit[j]=(long)p2;      for (j=1; j < ru; j++)
     for (i=1; i<=lastbid; i++) p2[i]=p1[i];        if (gcmp1(gcoeff(u,ru-1,j))) break;
     for (   ; i<=nblig; i++) p2[i]=zero;      if (j < ru)
     for (k=1; k<=lSml2; k++)      {
       mael(vecmunit2,k,j) = (long)zideallog(nf,(GEN)vunit[j],(GEN)vsml2[k]);        u = (GEN)u[j]; /* coords on (fu^ell, be) of a small generator */
     msign[j] = (long)zsigne(nf,(GEN)vunit[j],archartif);        ru--; setlg(u, ru);
         be = fix_be(bnfz, be, gmul(ell,u));
       }
   }    }
   for (j=1; j<=lSp; j++)    if (DEBUGLEVEL>1) fprintferr("beta LLL-reduced mod U^l = %Z\n",be);
   {    return reducebetanaive(bnfz, be, NULL, ell);
     p1 = zideallog(nf,(GEN)vecbeta[j],bid);  
     p2 = cgetg(nblig+1,t_COL); munit[j+nbvunit]=(long)p2;  
     for (i=1; i<=lastbid; i++) p2[i]=p1[i];  
     for (   ; i<=nblig; i++) p2[i]=coeff(matexpo,i-lastbid,j);  
     for (k=1; k<=lSml2; k++)  
       mael(vecmunit2,k,j+nbvunit) = (long)zideallog(nf,(GEN)vecbeta[j],(GEN)vsml2[k]);  
     msign[j+nbvunit] = (long)zsigne(nf,(GEN)vecbeta[j],archartif);  
   }  
   listx = get_listx(arch,msign,munit,vecmunit2,ell,lSp,nbvunit);  
   llistx= lg(listx);  
   listal= cgetg(llistx,t_VEC);  
   listg = concatsp(listg, concatsp(vecalpha0,vecbeta0));  
   l = lg(listg);  
   for (i=1; i<llistx; i++)  
   {  
     p1 = gun; p2 = (GEN)listx[i];  
     for (j=1; j<l; j++)  
       p1 = gmul(p1, powgi((GEN)listg[j],(GEN)p2[j]));  
     listal[i] = (long)reducealpha(nf,p1,gell);  
   }  
  /* A ce stade, tous les alpha dans la liste listal statisfont les  
   * congruences, noncongruences et signatures, et ne sont pas des puissances  
   * l-iemes donc x^l-alpha est irreductible, sa signature est correcte ainsi  
   * que son discriminant relatif. Reste a determiner son groupe de normes. */  
   if (DEBUGLEVEL) fprintferr("listalpha = %Z\n",listal);  
   p2 = cgetg(1,t_VEC);  
   for (i=1; i<llistx; i++)  
   {  
     p1 = gsub(gpuigs(polx[0],ell), (GEN)listal[i]);  
     if (all || gegal(rnfnormgroup(bnr,p1),subgroup)) p2 = concatsp(p2,p1);  
   }  
   if (all) return gcopy(p2);  
   switch(lg(p2)-1)  
   {  
     case 0: err(talker,"bug 6: no equation found in kummer");  
     case 1: break; /* OK */  
     default: fprintferr("equations = \n%Z",p2);  
       err(talker,"bug 7: more than one equation found in kummer");  
   }  
   return gcopy((GEN)p2[1]);  
 }  }
   
 static GEN  static GEN
 tauofideal(GEN id)  tauofalg(GEN x, GEN U) {
 {    return gsubst(lift(x), varn(U[1]), U);
   long j;  
   GEN p1,p2;  
   
   p1=gsubst(gmul((GEN)nfz[7],id),vnf,U);  
   p2=cgetg(lg(p1),t_MAT);  
   for (j=1; j<lg(p1); j++) p2[j]=(long)algtobasis(nfz,(GEN)p1[j]);  
   return p2;  
 }  }
   
 static GEN  static tau_s *
 tauofprimeideal(GEN pr)  get_tau(tau_s *tau, GEN nf, GEN U)
 {  {
   GEN p1 = dummycopy(pr);    GEN bas = (GEN)nf[7], Uzk;
     long i, l = lg(bas);
   p1[2] = (long)algtobasis(nfz, gsubst(gmul((GEN)nfz[7],(GEN)pr[2]),vnf,U));    Uzk = cgetg(l, t_MAT);
   return gcoeff(idealfactor(nfz,prime_to_ideal(nfz,p1)),1,1);    for (i=1; i<l; i++)
       Uzk[i] = (long)algtobasis(nf, tauofalg((GEN)bas[i], U));
     tau->zk = Uzk;
     tau->x  = U; return tau;
 }  }
   
 static long  static GEN tauoffamat(GEN x, tau_s *tau);
 isprimeidealconj(GEN pr1, GEN pr2)  
 {  
   GEN pr=pr1;  
   
   do  static GEN
   tauofelt(GEN x, tau_s *tau)
   {
     switch(typ(x))
   {    {
     if (gegal(pr,pr2)) return 1;      case t_COL: return gmul(tau->zk, x);
     pr = tauofprimeideal(pr);      case t_MAT: return tauoffamat(x, tau);
       default: return tauofalg(x, tau->x);
   }    }
   while (!gegal(pr,pr1));  
   return 0;  
 }  }
   static GEN
   tauofvec(GEN x, tau_s *tau)
   {
     long i, l = lg(x);
     GEN y = cgetg(l, typ(x));
   
 static long    for (i=1; i<l; i++) y[i] = (long)tauofelt((GEN)x[i], tau);
 isconjinprimelist(GEN listpr, GEN pr2)    return y;
   }
   /* [x, tau(x), ..., tau^m(x)] */
   static GEN
   powtau(GEN x, long m, tau_s *tau)
 {  {
   long ll=lg(listpr)-1,i;    GEN y = cgetg(m+1, t_VEC);
     long i;
     y[1] = (long)x;
     for (i=2; i<=m; i++) y[i] = (long)tauofelt((GEN)y[i-1], tau);
     return y;
   }
   
   for (i=1; i<=ll; i++)  static GEN
     if (isprimeidealconj((GEN)listpr[i],pr2)) return 1;  tauoffamat(GEN x, tau_s *tau)
   return 0;  {
     GEN y = cgetg(3, t_MAT);
     y[1] = (long)tauofvec((GEN)x[1], tau);
     y[2] = x[2]; return y;
 }  }
   
 static void  /* TODO: wasteful to reduce to 2-elt form. Compute image directly */
 computematexpoteta1(GEN A1, GEN R)  static GEN
   tauofideal(GEN nfz, GEN id, tau_s *tau)
 {  {
   long j;    GEN I = ideal_two_elt(nfz,id);
   GEN Aj = polun[vnf];    I[2] = (long)tauofelt((GEN)I[2], tau);
   matexpoteta1 = cgetg(degK+1,t_MAT);    return prime_to_ideal(nfz,I);
   for (j=1; j<=degK; j++)  }
   
   static int
   isprimeidealconj(GEN nfz, GEN pr1, GEN pr2, tau_s *tau)
   {
     GEN p = (GEN)pr1[1];
     GEN x = (GEN)pr1[2];
     GEN b1= (GEN)pr1[5];
     GEN b2= (GEN)pr2[5];
     if (!egalii(p, (GEN)pr2[1])
      || !egalii((GEN)pr1[3], (GEN)pr2[3])
      || !egalii((GEN)pr1[4], (GEN)pr2[4])) return 0;
     if (gegal(x,(GEN)pr2[2])) return 1;
     for(;;)
   {    {
     matexpoteta1[j] = (long)pol_to_vec(Aj, degKz);      if (int_elt_val(nfz,x,p,b2,NULL)) return 1;
     if (j<degK) Aj = gmod(gmul(Aj,A1), R);      x = FpV_red(tauofelt(x, tau), p);
       if (int_elt_val(nfz,x,p,b1,NULL)) return 0;
   }    }
 }  }
   
 static GEN  static int
 downtoK(GEN x)  isconjinprimelist(GEN nfz, GEN S, GEN pr, tau_s *tau)
 {  {
   long i;    long i, l;
   GEN p2,p3;  
     if (!tau) return 0;
     l = lg(S);
     for (i=1; i<l; i++)
       if (isprimeidealconj(nfz, (GEN)S[i],pr,tau)) return 1;
     return 0;
   }
   
   p2 = inverseimage(matexpoteta1, pol_to_vec(lift(x), degKz));  static GEN
   if (lg(p2)==1) err(talker,"not an element of K in downtoK");  downtoK(toK_s *T, GEN x)
   p3 = (GEN)p2[degK];  {
   for (i=degK-1; i; i--) p3 = gadd((GEN)p2[i],gmul(polx[vnf],p3));    long degKz = lg(T->invexpoteta1) - 1;
   return gmodulcp(p3,polnf);    GEN y = gmul(T->invexpoteta1, pol_to_vec(lift(x), degKz));
     return gmodulcp(gtopolyrev(y,varn(T->polnf)), T->polnf);
 }  }
   
 static GEN  static GEN
 tracetoK(GEN x)  tracetoK(toK_s *T, GEN x)
 {  {
     GEN p1 = x;
   long i;    long i;
   GEN p1,p2;    for (i=1; i < T->m; i++) p1 = gadd(x, tauofelt(p1,T->tau));
     return downtoK(T, p1);
   p1=x; p2=x;  
   for (i=1; i<=m-1; i++)  
   {  
     p1=gsubst(lift(p1),vnf,U);  
     p2=gadd(p2,p1);  
   }  
   return downtoK(p2);  
 }  }
   
 static GEN  static GEN
 normtoK(GEN x)  normtoK(toK_s *T, GEN x)
 {  {
     GEN p1 = x;
   long i;    long i;
   GEN p1,p2;    for (i=1; i < T->m; i++) p1 = gmul(x, tauofelt(p1,T->tau));
     return downtoK(T, p1);
   p1=x; p2=x;  
   for (i=1; i<=m-1; i++)  
   {  
     p1=gsubst(lift(p1),vnf,U);  
     p2=gmul(p2,p1);  
   }  
   return downtoK(p2);  
 }  }
   
 static GEN  static GEN
 computepolrel(void)  no_sol(long all, int i)
 {  {
   long i,j;    if (!all) err(talker,"bug%d in kummer",i);
   GEN p1,p2;    return cgetg(1,t_VEC);
   
   p1=gun; p2=gmodulcp(polx[vnf],R);  
   for (i=0; i<=m-1; i++)  
   {  
     p1=gmul(p1,gsub(polx[0],p2));  
     if (i<m-1) p2=gsubst(lift(p2),vnf,U);  
   }  
   for (j=2; j<=m+2; j++) p1[j]=(long)downtoK((GEN)p1[j]);  
   return p1;  
 }  }
   
 /* alg. 5.2.15. with remark */  
 static GEN  static GEN
 isprincipalell(GEN id)  get_gell(GEN bnr, GEN subgp, long all)
 {  {
   long i, l = lg(vecalpha);    GEN gell;
   GEN y,logdisc,be;    if (all)        gell = stoi(all);
     else if (subgp) gell = det(subgp);
   y = isprincipalgenforce(bnfz,id);    else            gell = det(diagonal(gmael(bnr,5,2)));
   logdisc = (GEN)y[1];    if (typ(gell) != t_INT) err(arither1);
   be = basistoalg(bnfz,(GEN)y[2]);    return gell;
   for (i=rc+1; i<l; i++)  
     be = gmul(be, powgi((GEN)vecalpha[i], modii(mulii((GEN)logdisc[i],(GEN)uu[i]),gell)));  
   y = cgetg(3,t_VEC);  
   y[1]=(long)logdisc; setlg(logdisc,rc+1);  
   y[2]=(long)be;  
   return y;  
 }  }
   
 /* alg. 5.3.11. */  typedef struct {
 static GEN    GEN Sm, Sml1, Sml2, Sl, ESml2;
 isvirtualunit(GEN v)  } primlist;
   
   static int
   build_list_Hecke(primlist *L, GEN nfz, GEN fa, GEN gothf, GEN gell, tau_s *tau)
 {  {
   long llist,i,j,ex, l = lg(vecalpha);    GEN listpr, listex, pr, p, factell;
   GEN p1,listex,listpr,q,ga,eps,vecy,logdisc;    long vd, vp, e, i, l, ell = itos(gell), degKz = degpol(nfz[1]);
   
   p1=idealfactor(nfz,v);    if (!fa) fa = idealfactor(nfz, gothf);
   listpr = (GEN)p1[1]; llist = lg(listpr)-1;    listpr = (GEN)fa[1];
   listex = (GEN)p1[2]; q = idmat(degKz);    listex = (GEN)fa[2]; l = lg(listpr);
   for (i=1; i<=llist; i++)    L->Sm  = cget1(l,t_VEC);
     L->Sml1= cget1(l,t_VEC);
     L->Sml2= cget1(l,t_VEC);
     L->Sl  = cget1(l+degKz,t_VEC);
     L->ESml2=cget1(l,t_VECSMALL);
     for (i=1; i<l; i++)
   {    {
     ex = itos((GEN)listex[i]);      pr = (GEN)listpr[i]; p = (GEN)pr[1]; e = itos((GEN)pr[3]);
     if (ex%ell) err(talker,"not a virtual unit in isvirtualunit");      vp = itos((GEN)listex[i]);
     q = idealmul(nfz,q, idealpows(nfz,(GEN)listpr[i], ex/ell));      if (!egalii(p,gell))
       {
         if (vp != 1) return 1;
         if (!isconjinprimelist(nfz, L->Sm,pr,tau)) appendL(L->Sm,pr);
       }
       else
       {
         vd = (vp-1)*(ell-1)-ell*e;
         if (vd > 0) return 4;
         if (vd==0)
         {
           if (!isconjinprimelist(nfz, L->Sml1,pr,tau)) appendL(L->Sml1, pr);
         }
         else
         {
           if (vp==1) return 2;
           if (!isconjinprimelist(nfz, L->Sml2,pr,tau))
           {
             appendL(L->Sml2, pr);
             appendL(L->ESml2,(GEN)vp);
           }
         }
       }
   }    }
   /* q^ell = (v) */    factell = primedec(nfz,gell); l = lg(factell);
   p1 = isprincipalgenforce(bnfz,q);    for (i=1; i<l; i++)
   logdisc = (GEN)p1[1];  
   ga = basistoalg(nfz,(GEN)p1[2]);  
   for (j=rc+1; j<l; j++)  
     ga = gmul(ga, powgi((GEN)vecalpha[j],divii((GEN)logdisc[j],(GEN)cyc[j])));  
   eps = gpuigs(ga,ell);  
   vecy = cgetg(rc+1,t_COL);  
   for (j=1; j<=rc; j++)  
   {    {
     vecy[j] = (long)divii((GEN)logdisc[j], divii((GEN)cyc[j],gell));      pr = (GEN)factell[i];
     eps = gmul(eps, powgi((GEN)vecalpha[j],(GEN)vecy[j]));      if (!idealval(nfz,gothf,pr))
         if (!isconjinprimelist(nfz, L->Sl,pr,tau)) appendL(L->Sl, pr);
   }    }
   eps = gdiv(v,eps);    return 0; /* OK */
   p1 = cgetg(3,t_VEC);  
   p1[1] = (long)concatsp(vecy, lift(isunit(bnfz,eps)));  
   p1[2] = (long)ga;  
   return p1;  
 }  }
   
 static GEN  static GEN
 lifttokz(GEN id, GEN A1)  logall(GEN nf, GEN vec, long lW, long mginv, long ell, GEN pr, long ex)
 {  {
   long i,j;    GEN m, M, bid = zidealstarinitgen(nf, idealpows(nf, pr, ex));
   GEN p1,p2,p3;    long ellrank, i, l = lg(vec);
   
   p1=gsubst(gmul((GEN)nf[7],id),vnf,A1);    ellrank = prank(gmael(bid,2,2), ell);
   p2=gmodulcp((GEN)nfz[7],R);    M = cgetg(l,t_MAT);
   p3=cgetg(degK*degKz+1,t_MAT);    for (i=1; i<l; i++)
   for (i=1; i<=degK; i++)  
     for (j=1; j<=degKz; j++)  
       p3[(i-1)*degKz+j]=(long)algtobasis(nfz,gmul((GEN)p1[i],(GEN)p2[j]));  
   return hnfmod(p3,detint(p3));  
 }  
   
 static GEN  
 steinitzaux(GEN id, GEN polrel)  
 {  
   long i,j;  
   GEN p1,p2,p3,vecid,matid,pseudomat,pid;  
   
   p1=gsubst(gmul((GEN)nfz[7],id),vnf,polx[0]);  
   for (j=1; j<=degKz; j++)  
     p1[j]=(long)gmod((GEN)p1[j],polrel);  
   p2=cgetg(degKz+1,t_MAT);  
   for (j=1; j<=degKz; j++)  
   {    {
     p3=cgetg(m+1,t_COL); p2[j]=(long)p3;      m = zideallog(nf, (GEN)vec[i], bid);
     for (i=1; i<=m; i++) p3[i]=(long)algtobasis(nf,truecoeff((GEN)p1[j],i-1));      setlg(m, ellrank+1);
       if (i < lW) m = gmulsg(mginv, m);
       M[i] = (long)m;
   }    }
   vecid=cgetg(degKz+1,t_VEC); matid=idmat(degK);    return M;
   for (j=1; j<=degKz; j++) vecid[j]=(long)matid;  
   pseudomat=cgetg(3,t_VEC);  
   pseudomat[1]=(long)p2; pseudomat[2]=(long)vecid;  
   pid=(GEN)nfhermite(nf,pseudomat)[2];  
   p1=matid;  
   for (j=1; j<=m; j++) p1=idealmul(nf,p1,(GEN)pid[j]);  
   return p1;  
 }  }
   
   /* compute the u_j (see remark 5.2.15.) */
 static GEN  static GEN
 normrelz(GEN id, GEN polrel, GEN steinitzZk)  get_u(GEN cyc, long rc, GEN gell)
 {  {
   GEN p1 = steinitzaux(idealhermite(nfz, id), polrel);    long i, l = lg(cyc);
   return idealdiv(nf,p1,steinitzZk);    GEN u = cgetg(l,t_VEC);
     for (i=1; i<=rc; i++) u[i] = zero;
     for (   ; i<  l; i++) u[i] = lmpinvmod((GEN)cyc[i], gell);
     return u;
 }  }
   
   /* alg. 5.2.15. with remark */
 static GEN  static GEN
 invimsubgroup(GEN bnrz, GEN bnr, GEN subgroup)  isprincipalell(GEN bnfz, GEN id, GEN cycgen, GEN u, GEN gell, long rc)
 {  {
   long l,j;    long i, l = lg(cycgen);
   GEN g,Plog,raycycz,rayclgpz,genraycycz,U,polrel,steinitzZk;    GEN logdisc, b, y = quick_isprincipalgen(bnfz, id);
   
   polrel = computepolrel();    logdisc = gmod((GEN)y[1], gell);
   steinitzZk = steinitzaux(idmat(degKz), polrel);    b = (GEN)y[2];
   rayclgpz = (GEN)bnrz[5];    for (i=rc+1; i<l; i++)
   raycycz   = (GEN)rayclgpz[2]; l=lg(raycycz);  
   genraycycz= (GEN)rayclgpz[3];  
   Plog = cgetg(l,t_MAT);  
   for (j=1; j<l; j++)  
   {    {
     g = normrelz((GEN)genraycycz[j],polrel,steinitzZk);      GEN e = modii(mulii((GEN)logdisc[i],(GEN)u[i]), gell);
     Plog[j] = (long)isprincipalray(bnr, g);      if (signe(e)) b = famat_mul(b, famat_pow((GEN)cycgen[i], e));
   }    }
   U = (GEN)hnfall(concatsp(Plog, subgroup))[2];    y = cgetg(3,t_VEC);
   setlg(U, l); for (j=1; j<l; j++) setlg(U[j], l);    y[1] = (long)logdisc; setlg(logdisc,rc+1);
   return hnfmod(concatsp(U, diagonal(raycycz)), (GEN)raycycz[1]);    y[2] = (long)b; return y;
 }  }
   
 static GEN  static GEN
 ideallogaux(long i, GEN al)  famat_factorback(GEN v, GEN e)
 {  {
   long llogli,valal;    long i, l = lg(e);
   GEN p1;    GEN V = cgetg(1, t_MAT);
     for (i=1; i<l; i++)
   valal = element_val(nfz,al,(GEN)listprSp[i]);      if (signe(e[i])) V = famat_mul(V, famat_pow((GEN)v[i], (GEN)e[i]));
   al = gmul(al,gpuigs((GEN)listunif[i],valal));    return V;
   p1 = zideallog(nfz,al,(GEN)listbid[i]);  
   llogli = listellrank[i];  
   setlg(p1,llogli+1); return p1;  
 }  }
   
 static GEN  static GEN
 ideallogauxsup(long i, GEN al)  compute_beta(GEN X, GEN vecWB, GEN ell, GEN bnfz)
 {  {
   long llogli,valal;    GEN BE, be;
   GEN p1;    BE = famat_reduce(famat_factorback(vecWB, X));
     BE[2] = (long)centermod((GEN)BE[2], ell);
     be = factorbackelt(BE, bnfz, NULL);
     be = reducebeta(bnfz, be, ell);
     if (DEBUGLEVEL>1) fprintferr("beta reduced = %Z\n",be);
     return basistoalg(bnfz, be); /* FIXME */
   }
   
   valal = element_val(nfz,al,(GEN)listprSp[i]);  static GEN
   al = gmul(al,gpuigs((GEN)listunif[i],valal));  get_Selmer(GEN bnf, GEN cycgen, long rc)
   p1 = zideallog(nfz,al,(GEN)listbidsup[i]);  {
   llogli = listellranksup[i];    GEN fu = check_units(bnf,"rnfkummer");
   setlg(p1,llogli+1); return p1;    GEN tu = gmael3(bnf,8,4,2);
     return concatsp(algtobasis(bnf,concatsp(fu,tu)), vecextract_i(cycgen,1,rc));
 }  }
   
   /* if all!=0, give all equations of degree 'all'. Assume bnr modulus is the
    * conductor */
 static GEN  static GEN
 vectau(GEN list)  rnfkummersimple(GEN bnr, GEN subgroup, GEN gell, long all)
 {  {
   long i,j,k,n;    long ell, i, j, degK, dK;
   GEN listz,listc,yz,yc,listfl,s, y = cgetg(3,t_VEC);    long lSml2, lSl2, lSp, rc, lW;
     long prec;
   
   listz = (GEN)list[1];    GEN bnf,nf,bid,ideal,arch,cycgen;
   listc = (GEN)list[2]; n = lg(listz);    GEN clgp,cyc;
   yz = cgetg(n,t_VEC); y[1] = (long)yz;    GEN Sp,listprSp,matP;
   yc = cgetg(n,t_VEC); y[2] = (long)yc;    GEN res,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
   listfl=cgetg(n,t_VECSMALL); for (i=1; i<n; i++) listfl[i] = 1;    primlist L;
   k = 1;  
   for (i=1; i<n; i++)    bnf = (GEN)bnr[1];
     nf  = (GEN)bnf[7];
     degK = degpol(nf[1]);
   
     bid = (GEN)bnr[2];
     ideal= gmael(bid,1,1);
     arch = gmael(bid,1,2); /* this is the conductor */
     ell = itos(gell);
     i = build_list_Hecke(&L, nf, (GEN)bid[3], ideal, gell, NULL);
     if (i) return no_sol(all,i);
   
     lSml2 = lg(L.Sml2)-1;
     Sp = concatsp(L.Sm, L.Sml1); lSp = lg(Sp)-1;
     listprSp = concatsp(L.Sml2, L.Sl); lSl2 = lg(listprSp)-1;
   
     cycgen = check_and_build_cycgen(bnf);
     clgp = gmael(bnf,8,1);
     cyc = (GEN)clgp[2]; rc = prank(cyc, ell);
   
     vecW = get_Selmer(bnf, cycgen, rc);
     u = get_u(cyc, rc, gell);
   
     vecBp = cgetg(lSp+1, t_VEC);
     matP  = cgetg(lSp+1, t_MAT);
     for (j=1; j<=lSp; j++)
   {    {
     if (!listfl[i]) continue;      GEN e, a, L;
       L = isprincipalell(bnf,(GEN)Sp[j], cycgen,u,gell,rc);
       e = (GEN)L[1]; matP[j]  = (long)e;
       a = (GEN)L[2]; vecBp[j] = (long)a;
     }
     vecWB = concatsp(vecW, vecBp);
   
     yz[k] = listz[i];    prec = DEFAULTPREC +
     s = (GEN)listc[i];      ((gexpo(vecWB) + gexpo(gmael(nf,5,1))) >> TWOPOTBYTES_IN_LONG);
     for (j=i+1; j<n; j++)    if (nfgetprec(nf) < prec) nf = nfnewprec(nf, prec);
     msign = zsigns(nf, vecWB);
   
     vecMsup = cgetg(lSml2+1,t_VEC);
     M = NULL;
     for (i=1; i<=lSl2; i++)
     {
       GEN pr = (GEN)listprSp[i];
       long e = itos((GEN)pr[3]), z = ell * (e / (ell-1));
   
       if (i <= lSml2)
     {      {
       if (listfl[j] && gegal((GEN)listz[j],(GEN)listz[i]))        z += 1 - L.ESml2[i];
       {        vecMsup[i] = (long)logall(nf, vecWB, 0,0, ell, pr,z+1);
         s = gadd(s,(GEN)listc[j]);  
         listfl[j] = 0;  
       }  
     }      }
     yc[k] = (long)s; k++;      M = vconcat(M, logall(nf, vecWB, 0,0, ell, pr,z));
   }    }
   setlg(yz, k);    lW = lg(vecW);
   setlg(yc, k); return y;    M = vconcat(M, concatsp(zeromat(rc,lW-1), matP));
 }  
   
 static GEN    K = FpM_ker(M, gell);
 subtau(GEN listx, GEN listy)    dK = lg(K)-1;
 {    y = cgetg(dK+1,t_VECSMALL);
   GEN y = cgetg(3,t_VEC);    res = cgetg(1,t_VEC); /* in case all = 1 */
   y[1] = (long)concatsp((GEN)listx[1], (GEN)listy[1]);    while (dK)
   y[2] = (long)concatsp((GEN)listx[2], gneg_i((GEN)listy[2]));    {
   return vectau(y);      for (i=1; i<dK; i++) y[i] = 0;
       y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
       do
       {
         GEN be, P, X = FpV_red(gmul_mati_smallvec(K, y), gell);
         if (ok_congruence(X, gell, lW, vecMsup) && ok_sign(X, msign, arch))
         {/* be satisfies all congruences, x^ell - be is irreducible, signature
           * and relative discriminant are correct */
           be = compute_beta(X, vecWB, gell, bnf);
           P = gsub(gpowgs(polx[0],ell), be);
           if (!all && gegal(rnfnormgroup(bnr,P),subgroup)) return P; /*DONE*/
           res = concatsp(res, P);
         }
       } while (increment(y, dK, ell));
       y[dK--] = 0;
     }
     if (all) return res;
     return gzero;
 }  }
   
   /* alg. 5.3.11 (return only discrete log mod ell) */
 static GEN  static GEN
 negtau(GEN list)  isvirtualunit(GEN bnf, GEN v, GEN cycgen, GEN cyc, GEN gell, long rc)
 {  {
   GEN y = cgetg(3,t_VEC);    GEN L, b, eps, y, q, nf = checknf(bnf);
   y[1] = list[1];    long i, l = lg(cycgen);
   y[2] = lneg((GEN)list[2]);  
   return y;    L = quick_isprincipalgen(bnf, idealsqrtn(nf, v, gell, 1));
     q = (GEN)L[1];
     if (gcmp0(q)) { eps = v; y = q; }
     else
     {
       b = (GEN)L[2];
       y = cgetg(l,t_COL);
       for (i=1; i<l; i++)
         y[i] = (long)diviiexact(mulii(gell,(GEN)q[i]), (GEN)cyc[i]);
       eps = famat_mul(famat_factorback(cycgen, y), famat_pow(b, gell));
       eps = famat_mul(famat_inv(eps), v);
     }
     setlg(y, rc+1);
     b = isunit(bnf,eps);
     if (lg(b) == 1) err(bugparier,"isvirtualunit");
     return concatsp(lift_intern(b), y);
 }  }
   
   /* id a vector of elements in nfz = relative extension of nf by polrel,
    * return the Steinitz element associated to the module generated by id */
 static GEN  static GEN
 multau(GEN listx, GEN listy)  steinitzaux(GEN nf, GEN id, GEN polrel)
 {  {
   GEN lzx,lzy,lcx,lcy,lzz,lcz, y = cgetg(3,t_VEC);    long i, degKz = lg(id)-1, degK = degpol(nf[1]);
   long nx,ny,i,j,k;    GEN V = dummycopy(id),vecid,matid,pseudomat,pid;
   
   lzx=(GEN)listx[1]; lcx=(GEN)listx[2]; nx=lg(lzx)-1;    vecid = cgetg(degKz+1,t_VEC); matid = idmat(degK);
   lzy=(GEN)listy[1]; lcy=(GEN)listy[2]; ny=lg(lzy)-1;    for (i=1; i<=degKz; i++) vecid[i] = (long)matid;
   lzz = cgetg(nx*ny+1,t_VEC); y[1]=(long)lzz;    for (i=1; i<=degKz; i++)
   lcz = cgetg(nx*ny+1,t_VEC); y[2]=(long)lcz;    {
   k = 0;      GEN v = (GEN)V[i];
   for (i=1; i<=nx; i++)      if (typ(v) == t_POL) { v = dummycopy(v); setvarn(v, 0); }
     for (j=1; j<=ny; j++)      V[i] = (long)gmod(v, polrel);
     {    }
       k++;    pseudomat = cgetg(3,t_VEC);
       lzz[k] = ladd((GEN)lzx[i],(GEN)lzy[j]);    pseudomat[1] = (long)vecpol_to_mat(V, degpol(polrel));
       lcz[k] = lmul((GEN)lcx[i],(GEN)lcy[j]);    pseudomat[2] = (long)vecid;
     }    pid = (GEN)nfhermite(nf,pseudomat)[2];
   return vectau(y);    return factorback(pid, nf); /* product */
 }  }
   
 static GEN  static GEN
 mulpoltau(GEN poltau, GEN list)  polrelKzK(toK_s *T, GEN x)
 {  {
   long i,j;    GEN P = roots_to_pol(powtau(x, T->m, T->tau), 0);
   GEN y;    long i, l = lg(P);
     for (i=2; i<l; i++) P[i] = (long)downtoK(T, (GEN)P[i]);
   j = lg(poltau)-2;    return P;
   y = cgetg(j+3,t_VEC);  
   y[1] = (long)negtau(multau(list,(GEN)poltau[1]));  
   for (i=2; i<=j+1; i++)  
     y[i] = (long)subtau((GEN)poltau[i-1],multau(list,(GEN)poltau[i]));  
   y[j+2] = poltau[j+1]; return y;  
 }  }
   
 /* th. 5.3.5. and prop. 5.3.9. */  /* N: Cl_m(Kz) --> Cl_m(K), lift subgroup from bnr to bnrz using Algo 4.1.11 */
 static GEN  static GEN
 computepolrelbeta(GEN be)  invimsubgroup(GEN bnrz, GEN bnr, GEN subgroup, toK_s *T)
 {  {
   long i,a,b,j;    long l, j;
   GEN e,u,u1,u2,u3,p1,p2,p3,p4,zet,be1,be2,listr,s,veczi,vecci,powtaubet;    GEN P,raycycz,rayclgpz,raygenz,U,polrel,steinitzZk;
     GEN nf = checknf(bnr), nfz = checknf(bnrz), polz = (GEN)nfz[1];
   
   switch (ell)    polrel = polrelKzK(T, polx[varn(polz)]);
     steinitzZk = steinitzaux(nf, (GEN)nfz[7], polrel);
     rayclgpz = (GEN)bnrz[5];
     raycycz = (GEN)rayclgpz[2]; l = lg(raycycz);
     raygenz = (GEN)rayclgpz[3];
     P = cgetg(l,t_MAT);
     for (j=1; j<l; j++)
   {    {
     case 2: err(talker,"you should not be here in rnfkummer !!"); break;      GEN g, id = idealhermite(nfz, (GEN)raygenz[j]);
     case 3: e=normtoK(be); u=tracetoK(be);      g = steinitzaux(nf, gmul((GEN)nfz[7], id), polrel);
       return gsub(gmul(polx[0],gsub(gsqr(polx[0]),gmulsg(3,e))),gmul(e,u));      g = idealdiv(nf, g, steinitzZk); /* N_{Kz/K}(gen[j]) */
     case 5: e=normtoK(be);      P[j] = (long)isprincipalray(bnr, g);
       if (d==2)  
       {  
         u=tracetoK(gpuigs(be,3));  
         p1=gadd(gmulsg(5,gsqr(e)), gmul(gsqr(polx[0]), gsub(gsqr(polx[0]),gmulsg(5,e))));  
         return gsub(gmul(polx[0],p1),gmul(e,u));  
       }  
       else  
       {  
         be1=gsubst(lift(be),vnf,U);  
         be2=gsubst(lift(be1),vnf,U);  
         u1=tracetoK(gmul(be,be1));  
         u2=tracetoK(gmul(gmul(be,be2),gsqr(be1)));  
         u3=tracetoK(gmul(gmul(gsqr(be),gpuigs(be1,3)),be2));  
         p1=gsub(gsqr(polx[0]),gmulsg(10,e));  
         p1=gsub(gmul(polx[0],p1),gmulsg(5,gmul(e,u1)));  
         p1=gadd(gmul(polx[0],p1),gmul(gmulsg(5,e),gsub(e,u2)));  
         p1=gsub(gmul(polx[0],p1),gmul(e,u3));  
         return p1;  
       }  
     default: p1=cgetg(2,t_VEC); p2=cgetg(3,t_VEC); p3=cgetg(2,t_VEC);  
       p4=cgetg(2,t_VEC); p3[1]=zero; p4[1]=un;  
       p2[1]=(long)p3; p2[2]=(long)p4; p1[1]=(long)p2;  
       zet=gmodulcp(polx[vnf], cyclo(ell,vnf));  
       listr=cgetg(m+1,t_VEC);  
       listr[1]=un;  
       for (i=2; i<=m; i++) listr[i]=(long)modii(mulii((GEN)listr[i-1],g),gell);  
       veczi=cgetg(m+1,t_VEC);  
       for (b=0; b<m; b++)  
       {  
         s=gzero;  
         for (a=0; a<m; a++)  
           s=gadd(gmul(polx[0],s),modii(mulii((GEN)listr[b+1],(GEN)listr[a+1]),gell));  
         veczi[b+1]=(long)s;  
       }  
       for (j=0; j<ell; j++)  
       {  
         vecci=cgetg(m+1,t_VEC);  
         for (b=0; b<m; b++) vecci[b+1]=(long)gpui(zet,mulsi(j,(GEN)listr[b+1]),0);  
         p4=cgetg(3,t_VEC);  
         p4[1]=(long)veczi; p4[2]=(long)vecci;  
         p1=mulpoltau(p1,p4);  
       }  
       powtaubet=cgetg(m+1,t_VEC);  
       powtaubet[1]=(long)be;  
       for (i=2; i<=m; i++) powtaubet[i]=(long)gsubst(lift((GEN)powtaubet[i-1]),vnf,U);  
       err(impl,"difficult Kummer for ell>=7"); return gzero;  
   }    }
   return NULL; /* not reached */    U = (GEN)hnfall(concatsp(P, subgroup))[2];
     setlg(U, l); for (j=1; j<l; j++) setlg(U[j], l);
     return hnfmod(concatsp(U, diagonal(raycycz)), (GEN)raycycz[1]);
 }  }
   
   /* t(b1,...,b_{k-1}) [prop 5.3.9] */
 static GEN  static GEN
 fix_be(GEN be, GEN u)  compute_t(GEN b, GEN r, long m, long ell)
 {  {
   GEN e,g, nf = checknf(bnfz), fu = gmael(bnfz,8,5);    long z, s, a, i, k = lg(b);
   long i,lu;    GEN t = cgetg(m+1, t_COL);
   
   lu = lg(u);    for (a = 0; a < m; a++)
   for (i=1; i<lu; i++)  
   {    {
     if (!signe(u[i])) continue;      z = m-1 - a; s = r[1 + z] - ell;
     e = mulsi(ell, (GEN)u[i]);      for (i=1; i < k; i++) s += r[1 + ((z + b[i]) % m)];
     g = gmodulcp((GEN)fu[i],(GEN)nf[1]);      t[1 + a] = lstoi(s / ell);
     be = gmul(be, powgi(g, e));  
   }    }
   return be;    return t;
 }  }
   
   /* Return multinomial(k-1; m1,...,m_{k-1}) where the mi are the
    * multiplicities of the bi [ assume b1 <= ... <= b_{k-1} ] */
 static GEN  static GEN
 logarch2arch(GEN x, long r1, long prec)  get_multinomial(GEN b)
 {  {
   long i, lx = lg(x), tx = typ(x);    long i, k = lg(b), a, m;
   GEN y = cgetg(lx, tx);    GEN z = mpfact(k - 1);
   
   if (tx == t_MAT)    a = b[1]; m = 1;
     for (i = 2; i < k; i++)
   {    {
     for (i=1; i<lx; i++) y[i] = (long)logarch2arch((GEN)x[i], r1, prec);      if (b[i] == a) m++;
     return y;      else
       {
         if (m > 1) z = diviiexact(z, mpfact(m));
         a = b[i]; m = 1;
       }
   }    }
   for (i=1; i<=r1;i++) y[i] = lexp((GEN)x[i],prec);    if (m > 1) z = diviiexact(z, mpfact(m));
   for (   ; i<lx; i++) y[i] = lexp(gmul2n((GEN)x[i],-1),prec);    return z;
   return y;  
 }  }
   
 /* multiply be by ell-th powers of units as to find small L2-norm for new be */  /* r[b[1]] + ... + r[b[k-1]] + 1 = 0 mod ell ?*/
   static int
   b_suitable(GEN b, GEN r, long k, long ell)
   {
     long i, s = 1;
     for (i = 1; i < k; i++) s += r[ 1 + b[i] ];
     return (s % ell) == 0;
   }
   
 static GEN  static GEN
 reducebetanaive(GEN be, GEN b)  pol_from_Newton(GEN S)
 {  {
   long i,k,n,ru,r1, prec = nfgetprec(bnfz);    long i, k, l = lg(S);
   GEN z,p1,p2,nmax,c, nf = checknf(bnfz);    GEN c = cgetg(l, t_VEC);
   
   if (DEBUGLEVEL) fprintferr("reduce modulo (Z_K^*)^l\n");    c[1] = S[1];
   r1 = nf_get_r1(nf);    for (k = 2; k < l; k++)
   if (!b) b = gmul(gmael(nf,5,1), algtobasis(nf,be));  
   n = max((ell>>1), 3);  
   z = cgetg(n+1, t_VEC);  
   c = gmulgs(greal((GEN)bnfz[3]), ell);  
   c = logarch2arch(c, r1, prec); /* = embeddings of fu^ell */  
   c = gprec_w(gnorm(c), DEFAULTPREC);  
   b = gprec_w(gnorm(b), DEFAULTPREC); /* need little precision */  
   z[1] = (long)concatsp(c, vecinv(c));  
   for (k=2; k<=n; k++) z[k] = (long) vecmul((GEN)z[1], (GEN)z[k-1]);  
   nmax = T2_from_embed_norm(b, r1);  
   ru = lg(c)-1; c = zerovec(ru);  
   for(;;)  
   {    {
     GEN B = NULL;      GEN s = (GEN)S[k];
     long besti = 0, bestk = 0;      for (i = 1; i < k; i++) s = gadd(s, gmul((GEN)S[i], (GEN)c[k-i]));
     for (k=1; k<=n; k++)      c[k] = ldivgs(s, -k);
       for (i=1; i<=ru; i++)  
       {  
         p1 = vecmul(b, gmael(z,k,i));    p2 = T2_from_embed_norm(p1,r1);  
         if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk = k; continue; }  
         p1 = vecmul(b, gmael(z,k,i+ru)); p2 = T2_from_embed_norm(p1,r1);  
         if (gcmp(p2,nmax) < 0) { B=p1; nmax=p2; besti=i; bestk =-k; }  
       }  
     if (!B) break;  
     b = B; c[besti] = laddis((GEN)c[besti], bestk);  
   }    }
   if (DEBUGLEVEL) fprintferr("unit exponents = %Z\n",c);    return gadd(gpowgs(polx[0], l-1), gtopoly(c, 0));
   return fix_be(be,c);  
 }  }
   
   /* th. 5.3.5. and prop. 5.3.9. */
 static GEN  static GEN
 reducebeta(GEN be)  compute_polrel(GEN nfz, toK_s *T, GEN be, long g, long ell)
 {  {
   long j,ru, prec = nfgetprec(bnfz);    long i, k, m = T->m;
   GEN emb,z,u,matunit, nf = checknf(bnfz);    GEN r, powtaubet, S, X = polx[0], e = normtoK(T,be);
   
   if (gcmp1(gnorm(be))) return reducebetanaive(be,NULL);    switch (ell)
   matunit = gmulgs(greal((GEN)bnfz[3]), ell); /* log. embeddings of fu^ell */    { /* special-cased for efficiency */
   for (;;)      GEN p1, u;
       case 2: err(bugparier,"rnfkummer (-1 not in nf?)"); break;
       case 3: u = tracetoK(T,be);
         p1 = gsub(gsqr(X), gmulsg(3,e));
         return gsub(gmul(X,p1), gmul(e,u));
       case 5:
         if (m == 2)
         {
           u = tracetoK(T, gpowgs(be,3));
           p1 = gadd(gmulsg(5,gsqr(e)), gmul(gsqr(X), gsub(gsqr(X),gmulsg(5,e))));
           return gsub(gmul(X,p1), gmul(e,u));
         }
         else
         { /* m = 4 */
           GEN be1, be2, u1, u2, u3;
           be1 = tauofelt(be, T->tau);
           be2 = tauofelt(be1,T->tau);
           u1 = tracetoK(T, gmul(be,be1));
           u2 = tracetoK(T, gmul(gmul(be,be2),gsqr(be1)));
           u3 = tracetoK(T, gmul(gmul(gsqr(be),gpowgs(be1,3)),be2));
           p1 = gsub(gsqr(X), gmulsg(10,e));
           p1 = gsub(gmul(X,p1), gmulsg(5,gmul(e,u1)));
           p1 = gadd(gmul(X,p1), gmul(gmulsg(5,e),gsub(e,u2)));
           return gsub(gmul(X,p1), gmul(e,u3));
         }
     }
     /* general case */
     r = cgetg(m+1,t_VECSMALL); /* r[i+1] = g^i mod ell */
     r[1] = 1;
     for (i=2; i<=m; i++) r[i] = (r[i-1] * g) % ell;
     powtaubet = powtau(be, m, T->tau);
     S = cgetg(ell+1, t_VEC); /* Newton sums */
     S[1] = zero;
     for (k = 2; k <= ell; k++)
   {    {
     z = get_arch_real(nf, be, &emb, prec);      GEN z, g = gzero, b = vecsmall_const(k-1, 0);
     if (z) break;      do
     prec = (prec-1)<<1;      {
     if (DEBUGLEVEL) err(warnprec,"reducebeta",prec);        if (! b_suitable(b, r, k, ell)) continue;
     nf = nfnewprec(nf,prec);        z = factorbackelt(powtaubet, compute_t(b, r, m, ell), nfz);
         if (typ(z) == t_COL) z = basistoalg(nfz, z);
         g = gadd(g, gmul(get_multinomial(b), z));
       } while (increment_inc(b, k, m));
       S[k] = lmul(gmulsg(ell, e), tracetoK(T,g));
   }    }
   z = concatsp(matunit, z);    return pol_from_Newton(S);
   u = lllintern(z, 1, prec);  
   if (!u) return reducebetanaive(be,emb); /* shouldn't occur */  
   ru = lg(u);  
   for (j=1; j<ru; j++)  
     if (smodis(gcoeff(u,ru-1,j), ell)) break; /* prime to ell */  
   u = (GEN)u[j]; /* coords on (fu^ell, be) of a small generator */  
   ru--; setlg(u, ru);  
   be = powgi(be, (GEN)u[ru]);  
   return reducebetanaive(fix_be(be, u), NULL);  
 }  }
   
 /* cf. algo 5.3.18 */  typedef struct {
     GEN R; /* compositum(P,Q) */
     GEN p; /* Mod(p,R) root of P */
     GEN q; /* Mod(q,R) root of Q */
     GEN k; /* Q[X]/R generated by q + k p */
     GEN rev;
   } compo_s;
   
 static GEN  static GEN
 testx(GEN bnf, GEN X, GEN module, GEN subgroup, GEN vecMsup)  lifttoKz(GEN nfz, GEN nf, GEN id, compo_s *C)
 {  {
   long i,v,l,lX;    GEN I = ideal_two_elt(nf,id);
   GEN be,polrelbe,p1;    GEN x = gmul((GEN)nf[7], (GEN)I[2]);
     I[2] = (long)algtobasis(nfz, RX_RXQ_compo(x, C->p, C->R));
   if (gcmp0(X)) return NULL;    return prime_to_ideal(nfz,I);
   lX = lg(X);  
   for (i=dv+1; i<lX; i++)  
     if (gcmp0((GEN)X[i])) return NULL;  
   l = lg(vecMsup);  
   for (i=1; i<l; i++)  
     if (gcmp0(FpV_red(gmul((GEN)vecMsup[i],X), gell))) return NULL;  
   be = gun;  
   for (i=1; i<lX; i++)  
     be = gmul(be, powgi((GEN)vecw[i], (GEN)X[i]));  
   if (DEBUGLEVEL>1) fprintferr("reducing beta = %Z\n",be);  
   be = reducebeta(be);  
   if (DEBUGLEVEL>1) fprintferr("beta reduced = %Z\n",be);  
   polrelbe = computepolrelbeta(be);  
   v = varn(polrelbe);  
   p1 = unifpol((GEN)bnf[7],polrelbe,0);  
   p1 = denom(gtovec(p1));  
   polrelbe = gsubst(polrelbe,v, gdiv(polx[v],p1));  
   polrelbe = gmul(polrelbe, gpowgs(p1, degpol(polrelbe)));  
   if (DEBUGLEVEL>1) fprintferr("polrelbe = %Z\n",polrelbe);  
   p1 = rnfconductor(bnf,polrelbe,0);  
   if (!gegal((GEN)p1[1],module) || !gegal((GEN)p1[3],subgroup)) return NULL;  
   return polrelbe;  
 }  }
   
   static void
   compositum_red(compo_s *C, GEN P, GEN Q)
   {
     GEN a, z = (GEN)compositum2(P, Q)[1];
     C->R = (GEN)z[1];
     C->p = (GEN)z[2];
     C->q = (GEN)z[3];
     C->k = (GEN)z[4];
     /* reduce R */
     z = polredabs0(C->R, nf_ORIG|nf_PARTIALFACT);
     C->R = (GEN)z[1];
     if (DEBUGLEVEL>1) fprintferr("polred(compositum) = %Z\n",C->R);
     a    = (GEN)z[2];
     C->p = poleval(lift_intern(C->p), a);
     C->q = poleval(lift_intern(C->q), a);
     C->rev = modreverse_i((GEN)a[2], (GEN)a[1]);
   }
   
 GEN  static GEN
 rnfkummer(GEN bnr, GEN subgroup, long all, long prec)  _rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
 {  {
   long i,j,l,av=avma,e,vp,vd,dK,lSl,lSp,lSl2,lSml2,dc,ru,rv,nbcol;    long ell, i, j, m, d, dK, dc, rc, ru, rv, g, mginv, degK, degKz, vnf;
   GEN p1,p2,p3,p4,wk,pr;    long l, lSp, lSml2, lSl2, lW;
   GEN bnf,rayclgp,bid,ideal,cycgen,vselmer;    GEN polnf,bnf,nf,bnfz,nfz,bid,ideal,cycgen,gell,p1,p2,wk,U,vselmer;
   GEN kk,clgp,fununits,torsunit,vecB,vecC,Tc,Tv,P;    GEN clgp,cyc,gen;
   GEN Q,Qt,idealz,gothf,factgothf,listpr,listex,factell,p,vecnul;    GEN Q,idealz,gothf;
   GEN M,al,K,Msup,X,finalresult,y,module,A1,A2,vecMsup;    GEN res,u,M,K,y,vecMsup,vecW,vecWA,vecWB,vecB,vecC,vecAp,vecBp;
   GEN listmodsup,vecalphap,vecbetap,mginv,matP,ESml2,Sp,Sm,Sml1,Sml2,Sl;    GEN matP,Sp,listprSp,Tc,Tv,P;
     primlist L;
     toK_s T;
     tau_s _tau, *tau;
     compo_s COMPO;
   
   checkbnrgen(bnr);    checkbnrgen(bnr);
   wk = gmael4(bnr,1,8,4,1);  
   if (all) gell = stoi(all);  
   else  
   {  
     if (!gcmp0(subgroup)) gell = det(subgroup);  
     else gell = det(diagonal(gmael(bnr,5,2)));  
   }  
   if (gcmp1(gell)) { avma = av; return polx[varn(gmael3(bnr,1,7,1))]; }  
   if (!isprime(gell)) err(impl,"kummer for composite relative degree");  
   if (divise(wk,gell))  
     return gerepileupto(av,rnfkummersimple(bnr,subgroup,all));  
   if (all && gcmp0(subgroup))  
     err(talker,"kummer when zeta not in K requires a specific subgroup");  
   bnf = (GEN)bnr[1];    bnf = (GEN)bnr[1];
   nf = (GEN)bnf[7];    nf  = (GEN)bnf[7];
   polnf = (GEN)nf[1]; vnf = varn(polnf); degK = degpol(polnf);    polnf = (GEN)nf[1]; vnf = varn(polnf);
   if (!vnf) err(talker,"main variable in kummer must not be x");    if (!vnf) err(talker,"main variable in kummer must not be x");
       /* step 7 */    wk = gmael3(bnf,8,4,1);
   p1 = conductor(bnr,subgroup,2);    /* step 7 */
       /* fin step 7 */    if (all) subgroup = NULL;
   bnr = (GEN)p1[2];    p1 = conductor(bnr, subgroup, 2);
     bnr      = (GEN)p1[2];
   subgroup = (GEN)p1[3];    subgroup = (GEN)p1[3];
   rayclgp = (GEN)bnr[5];    gell = get_gell(bnr,subgroup,all);
   raycyc = (GEN)rayclgp[2];    if (gcmp1(gell)) return polx[vnf];
   bid = (GEN)bnr[2]; module=(GEN)bid[1];    if (!isprime(gell)) err(impl,"kummer for composite relative degree");
   ideal = (GEN)module[1];    if (divise(wk,gell)) return rnfkummersimple(bnr, subgroup, gell, all);
   if (gcmp0(subgroup)) subgroup = diagonal(raycyc);  
     bid = (GEN)bnr[2];
     ideal = gmael(bid,1,1);
   ell = itos(gell);    ell = itos(gell);
       /* step 1 of alg 5.3.5. */    /* step 1 of alg 5.3.5. */
   if (DEBUGLEVEL>2) fprintferr("Step 1\n");    if (DEBUGLEVEL>2) fprintferr("Step 1\n");
     compositum_red(&COMPO, polnf, cyclo(ell,vnf));
   p1 = (GEN)compositum2(polnf, cyclo(ell,vnf))[1];    /* step 2 */
   R = (GEN)p1[1];  
   A1= (GEN)p1[2];  
   A2= (GEN)p1[3];  
   kk= (GEN)p1[4];  
   if (signe(leadingcoeff(R)) < 0)  
   {  
     R = gneg_i(R);  
     A1 = gmodulcp(lift(A1),R);  
     A2 = gmodulcp(lift(A2),R);  
   }  
       /* step 2 */  
   if (DEBUGLEVEL>2) fprintferr("Step 2\n");    if (DEBUGLEVEL>2) fprintferr("Step 2\n");
   degKz = degpol(R);    degK  = degpol(polnf);
   m = degKz/degK;    degKz = degpol(COMPO.R);
   d = (ell-1)/m;    m = degKz / degK;
   g = lift(gpuigs(gener(gell),d)); /* g has order m in all (Z/ell^k)^* */    d = (ell-1) / m;
   if (gcmp1(powmodulo(g, stoi(m), stoi(ell*ell)))) g = addsi(ell,g);    g = powuumod(u_gener(ell), d, ell);
       /* step reduction of R */    if (powuumod(g, m, ell*ell) == 1) g += ell; /* ord(g)=m in all (Z/ell^k)^* */
   if (degKz<=20)    /* step 3 */
   {  
     GEN A3,A3rev;  
   
     if (DEBUGLEVEL>2) fprintferr("Step reduction\n");  
     p1 = polredabs2(R,prec);  
     if (DEBUGLEVEL>2) fprintferr("polredabs = %Z",p1[1]);  
     R = (GEN)p1[1];  
     A3= (GEN)p1[2];  
     A1 = poleval(lift(A1), A3);  
     A2 = poleval(lift(A2), A3);  
     A3rev= polymodrecip(A3);  
     U = gadd(powgi(A2,g), gmul(kk,A1));  
     U = poleval(lift(A3rev), U);  
   }  
   else U = gadd(powgi(A2,g), gmul(kk,A1));  
       /* step 3 */  
       /* one could factor disc(R) using th. 2.1.6. */  
   if (DEBUGLEVEL>2) fprintferr("Step 3\n");    if (DEBUGLEVEL>2) fprintferr("Step 3\n");
   bnfz = bnfinit0(R,1,NULL,prec);    /* could factor disc(R) using th. 2.1.6. */
     bnfz = bnfinit0(COMPO.R,1,NULL,prec);
     cycgen = check_and_build_cycgen(bnfz);
   nfz = (GEN)bnfz[7];    nfz = (GEN)bnfz[7];
   clgp = gmael(bnfz,8,1);    clgp = gmael(bnfz,8,1);
   cyc   = (GEN)clgp[2]; rc = ellrank(cyc,ell);    cyc = (GEN)clgp[2]; rc = prank(cyc,ell);
   gencyc= (GEN)clgp[3]; l = lg(cyc);    gen = (GEN)clgp[3]; l = lg(cyc);
   vecalpha=cgetg(l,t_VEC);    u = get_u(cyc, rc, gell);
   cycgen = check_and_build_cycgen(bnfz);  
   for (j=1; j<l; j++)  
     vecalpha[j] = (long)basistoalg(nfz, famat_to_nf(nfz, (GEN)cycgen[j]));  
       /* computation of the uu(j) (see remark 5.2.15.) */  
   uu = cgetg(l,t_VEC);  
   for (j=1; j<=rc; j++) uu[j] = zero;  
   for (   ; j<  l; j++) uu[j] = lmpinvmod((GEN)cyc[j], gell);  
   
   fununits = check_units(bnfz,"rnfkummer");    vselmer = get_Selmer(bnfz, cycgen, rc);
   torsunit = gmael3(bnfz,8,4,2);    ru = (degKz>>1)-1;
   ru=(degKz>>1)-1;    rv = rc+ru+1;
   rv=rc+ru+1;  
   vselmer = cgetg(rv+1,t_VEC);    /* compute action of tau */
   for (j=1; j<=rc; j++) vselmer[j] = vecalpha[j];    U = gadd(gpowgs(COMPO.q, g), gmul(COMPO.k, COMPO.p));
   for (   ; j< rv; j++) vselmer[j] = lmodulcp((GEN)fununits[j-rc],R);    U = poleval(COMPO.rev, U);
   vselmer[rv]=(long)gmodulcp((GEN)torsunit,R);    tau = get_tau(&_tau, nfz, U);
     /* step 4 */  
     /* step 4 */
   if (DEBUGLEVEL>2) fprintferr("Step 4\n");    if (DEBUGLEVEL>2) fprintferr("Step 4\n");
   vecB=cgetg(rc+1,t_VEC);    vecB=cgetg(rc+1,t_VEC);
   Tc=cgetg(rc+1,t_MAT);    Tc=cgetg(rc+1,t_MAT);
   for (j=1; j<=rc; j++)    for (j=1; j<=rc; j++)
   {    {
     p1 = isprincipalell(tauofideal((GEN)gencyc[j]));      p1 = tauofideal(nfz, (GEN)gen[j], tau);
       p1 = isprincipalell(bnfz, p1, cycgen,u,gell,rc);
     Tc[j]  = p1[1];      Tc[j]  = p1[1];
     vecB[j]= p1[2];      vecB[j]= p1[2];
   }    }
   p1=cgetg(m,t_VEC);  
   p1[1]=(long)idmat(rc);    vecC = cgetg(rc+1,t_VEC);
   for (j=2; j<=m-1; j++) p1[j]=lmul((GEN)p1[j-1],Tc);    for (j=1; j<=rc; j++) vecC[j] = lgetg(1, t_MAT);
   p2=cgetg(rc+1,t_VEC);    p1 = cgetg(m,t_VEC);
   for (j=1; j<=rc; j++) p2[j]=un;    p1[1] = (long)idmat(rc);
   p3=vecB;    for (j=2; j<=m-1; j++) p1[j] = lmul((GEN)p1[j-1],Tc);
     p2 = vecB;
   for (j=1; j<=m-1; j++)    for (j=1; j<=m-1; j++)
   {    {
     p3 = gsubst(lift(p3),vnf,U);      GEN T = FpM_red(gmulsg((j*d)%ell,(GEN)p1[m-j]), gell);
     p4 = groupproduct(grouppows(p3,(j*d)%ell),(GEN)p1[m-j]);      p2 = tauofvec(p2, tau);
     for (i=1; i<=rc; i++) p2[i] = lmul((GEN)p2[i],(GEN)p4[i]);      for (i=1; i<=rc; i++)
         vecC[i] = (long)famat_mul((GEN)vecC[i], famat_factorback(p2, (GEN)T[i]));
   }    }
   vecC=p2;    for (i=1; i<=rc; i++) vecC[i] = (long)famat_reduce((GEN)vecC[i]);
       /* step 5 */    /* step 5 */
   if (DEBUGLEVEL>2) fprintferr("Step 5\n");    if (DEBUGLEVEL>2) fprintferr("Step 5\n");
   Tv = cgetg(rv+1,t_MAT);    Tv = cgetg(rv+1,t_MAT);
   for (j=1; j<=rv; j++)    for (j=1; j<=rv; j++)
   {    {
     Tv[j] = isvirtualunit(gsubst(lift((GEN)vselmer[j]),vnf,U))[1];      p1 = tauofelt((GEN)vselmer[j], tau);
     if (DEBUGLEVEL>2) fprintferr("   %ld\n",j);      if (typ(p1) == t_MAT) /* famat */
         p1 = factorbackelt((GEN)p1[1], FpV_red((GEN)p1[2],gell), nfz);
       Tv[j] = (long)isvirtualunit(bnfz, p1, cycgen,cyc,gell,rc);
   }    }
   P = FpM_ker(gsub(Tv, g), gell);    P = FpM_ker(gsubgs(Tv, g), gell);
   dv= lg(P)-1; vecw = cgetg(dv+1,t_VEC);    lW = lg(P); vecW = cgetg(lW,t_VEC);
   for (j=1; j<=dv; j++)    for (j=1; j<lW; j++) vecW[j] = (long)famat_factorback(vselmer, (GEN)P[j]);
   {    /* step 6 */
     p1 = gun;  
     for (i=1; i<=rv; i++) p1 = gmul(p1, powgi((GEN)vselmer[i],gcoeff(P,i,j)));  
     vecw[j] = (long)p1;  
   }  
       /* step 6 */  
   if (DEBUGLEVEL>2) fprintferr("Step 6\n");    if (DEBUGLEVEL>2) fprintferr("Step 6\n");
   Q = FpM_ker(gsub(gtrans(Tc), g), gell);    Q = FpM_ker(gsubgs(gtrans(Tc), g), gell);
   Qt = gtrans(Q); dc = lg(Q)-1;    /* step 8 */
       /* step 7 done above */    if (DEBUGLEVEL>2) fprintferr("Step 8\n");
       /* step 8 */    p1 = RXQ_powers(lift_intern(COMPO.p), COMPO.R, degK-1);
   if (DEBUGLEVEL>2) fprintferr("Step 7 and 8\n");    p1 = vecpol_to_mat(p1, degKz);
   idealz=lifttokz(ideal, A1);    T.invexpoteta1 = invmat(p1); /* left inverse */
   computematexpoteta1(lift(A1), R);    T.polnf = polnf;
   if (!divise(idealnorm(nf,ideal),gell)) gothf = idealz;    T.tau = tau;
     T.m = m;
   
     idealz = lifttoKz(nfz, nf, ideal, &COMPO);
     if (smodis(idealnorm(nf,ideal), ell)) gothf = idealz;
   else    else
   {    { /* ell | N(ideal) */
     GEN bnrz, subgroupz;      GEN bnrz = buchrayinitgen(bnfz, idealz);
     bnrz = buchrayinitgen(bnfz,idealz);      GEN subgroupz = invimsubgroup(bnrz, bnr, subgroup, &T);
     subgroupz = invimsubgroup(bnrz,bnr,subgroup);  
     gothf = conductor(bnrz,subgroupz,0);      gothf = conductor(bnrz,subgroupz,0);
   }    }
       /* step 9 */    /* step 9, 10, 11 */
   if (DEBUGLEVEL>2) fprintferr("Step 9\n");    if (DEBUGLEVEL>2) fprintferr("Step 9, 10 and 11\n");
   factgothf=idealfactor(nfz,gothf);    i = build_list_Hecke(&L, nfz, NULL, gothf, gell, tau);
   listpr = (GEN)factgothf[1];    if (i) return no_sol(all,i);
   listex = (GEN)factgothf[2]; l = lg(listpr);  
       /* step 10 and step 11 */    lSml2 = lg(L.Sml2)-1;
   if (DEBUGLEVEL>2) fprintferr("Step 10 and 11\n");    Sp = concatsp(L.Sm, L.Sml1); lSp = lg(Sp)-1;
   Sm  = cgetg(l,t_VEC); setlg(Sm,1);    listprSp = concatsp(L.Sml2, L.Sl); lSl2 = lg(listprSp)-1;
   Sml1= cgetg(l,t_VEC); setlg(Sml1,1);  
   Sml2= cgetg(l,t_VEC); setlg(Sml2,1);    /* step 12 */
   Sl  = cgetg(l+degKz,t_VEC); setlg(Sl,1);  
   ESml2=cgetg(l,t_VECSMALL); setlg(ESml2,1);  
   for (i=1; i<l; i++)  
   {  
     pr = (GEN)listpr[i]; p = (GEN)pr[1]; e = itos((GEN)pr[3]);  
     vp = itos((GEN)listex[i]);  
     if (!egalii(p,gell))  
     {  
       if (vp != 1) { avma = av; return gzero; }  
       if (!isconjinprimelist(Sm,pr)) _append(Sm,pr);  
     }  
     else  
     {  
       vd = (vp-1)*(ell-1)-ell*e;  
       if (vd > 0) { avma = av; return gzero; }  
       if (vd==0)  
       {  
         if (!isconjinprimelist(Sml1,pr)) _append(Sml1, pr);  
       }  
       else  
       {  
         if (vp==1) { avma = av; return gzero; }  
         if (!isconjinprimelist(Sml2,pr))  
         {  
           _append(Sml2, pr);  
           _append(ESml2,(GEN)vp);  
         }  
       }  
     }  
   }  
   factell = primedec(nfz,gell); l = lg(factell);  
   for (i=1; i<l; i++)  
   {  
     pr = (GEN)factell[i];  
     if (!idealval(nfz,gothf,pr))  
       if (!isconjinprimelist(Sl,pr)) _append(Sl, pr);  
   }  
   lSml2=lg(Sml2)-1; lSl=lg(Sl)-1; lSl2=lSml2+lSl;  
   Sp=concatsp(Sm,Sml1); lSp=lg(Sp)-1;  
       /* step 12 */  
   if (DEBUGLEVEL>2) fprintferr("Step 12\n");    if (DEBUGLEVEL>2) fprintferr("Step 12\n");
   vecbetap=cgetg(lSp+1,t_VEC);    vecAp = cgetg(lSp+1, t_VEC);
   vecalphap=cgetg(lSp+1,t_VEC);    vecBp = cgetg(lSp+1, t_VEC);
   matP=cgetg(lSp+1,t_MAT);    matP  = cgetg(lSp+1, t_MAT);
   for (j=1; j<=lSp; j++)    for (j=1; j<=lSp; j++)
   {    {
     p1=isprincipalell((GEN)Sp[j]);      GEN e, a, ap;
     matP[j]=p1[1];      p1 = isprincipalell(bnfz, (GEN)Sp[j], cycgen,u,gell,rc);
     p2=gun;      e = (GEN)p1[1]; matP[j] = (long)e;
     for (i=1; i<=rc; i++)      a = (GEN)p1[2];
       p2 = gmul(p2, powgi((GEN)vecC[i],gmael(p1,1,i)));      p2 = famat_mul(famat_factorback(vecC, gneg(e)), a);
     p3=gdiv((GEN)p1[2],p2); vecbetap[j]=(long)p3;      vecBp[j] = (long)p2;
     p2=gun;      ap = cgetg(1, t_MAT);
     for (i=0; i<m; i++)      for (i=0; i<m; i++)
     {      {
       p2 = gmul(p2, powgi(p3, powmodulo(g,stoi(m-1-i),gell)));        ap = famat_mul(ap, famat_pow(p2, utoi(powuumod(g,m-1-i,ell))));
       if (i<m-1) p3=gsubst(lift(p3),vnf,U);        if (i < m-1) p2 = tauofelt(p2, tau);
     }      }
     vecalphap[j]=(long)p2;      vecAp[j] = (long)ap;
   }    }
       /* step 13 */    /* step 13 */
   if (DEBUGLEVEL>2) fprintferr("Step 13\n");    if (DEBUGLEVEL>2) fprintferr("Step 13\n");
   nbcol=lSp+dv;    vecWA = concatsp(vecW, vecAp);
   vecw=concatsp(vecw,vecbetap);    vecWB = concatsp(vecW, vecBp);
   listmod=cgetg(lSl2+1,t_VEC);  
   listunif=cgetg(lSl2+1,t_VEC);    /* step 14, 15, and 17 */
   listprSp = concatsp(Sml2, Sl);    if (DEBUGLEVEL>2) fprintferr("Step 14, 15 and 17\n");
   for (j=1; j<=lSl2; j++)    mginv = (m * u_invmod(g,ell)) % ell;
     vecMsup = cgetg(lSml2+1,t_VEC);
     M = NULL;
     for (i=1; i<=lSl2; i++)
   {    {
     GEN pr = (GEN)listprSp[j];      GEN pr = (GEN)listprSp[i];
     long e = itos((GEN)pr[3]), z = ell * (e / (ell-1));      long e = itos((GEN)pr[3]), z = ell * (e / (ell-1));
     if (j <= lSml2) z += 1 - ESml2[j];  
     listmod[j]=(long)idealpows(nfz,pr, z);      if (i <= lSml2)
     listunif[j]=(long)basistoalg(nfz,gdiv((GEN)pr[5],(GEN)pr[1]));      {
         z += 1 - L.ESml2[i];
         vecMsup[i] = (long)logall(nfz, vecWA,lW,mginv,ell, pr,z+1);
       }
       M = vconcat(M, logall(nfz, vecWA,lW,mginv,ell, pr,z));
   }    }
       /* step 14 and step 15 */    dc = lg(Q)-1;
   if (DEBUGLEVEL>2) fprintferr("Step 14 and 15\n");    if (dc)
   listbid=cgetg(lSl2+1,t_VEC);  
   listellrank = cgetg(lSl2+1,t_VECSMALL);  
   for (i=1; i<=lSl2; i++)  
   {    {
     listbid[i]=(long)zidealstarinitgen(nfz,(GEN)listmod[i]);      GEN QtP = gmul(gtrans_i(Q), matP);
     listellrank[i] = ellrank(gmael3(listbid,i,2,2), ell);      M = vconcat(M, concatsp(zeromat(dc,lW-1), QtP));
   }    }
   mginv = modii(mulsi(m, mpinvmod(g,gell)), gell);    if (!M) M = zeromat(1, lSp + lW - 1);
   vecnul=cgetg(dc+1,t_COL); for (i=1; i<=dc; i++) vecnul[i]=zero;    /* step 16 */
   M=cgetg(nbcol+1,t_MAT);  
   for (j=1; j<=dv; j++)  
   {  
     p1=cgetg(1,t_COL);  
     al=(GEN)vecw[j];  
     for (i=1; i<=lSl2; i++) p1 = concatsp(p1,ideallogaux(i,al));  
     p1=gmul(mginv,p1);  
     M[j]=(long)concatsp(p1,vecnul);  
   }  
   for (   ; j<=nbcol; j++)  
   {  
     p1=cgetg(1,t_COL);  
     al=(GEN)vecalphap[j-dv];  
     for (i=1; i<=lSl2; i++) p1 = concatsp(p1,ideallogaux(i,al));  
     M[j]=(long)concatsp(p1,gmul(Qt,(GEN)matP[j-dv]));  
   }  
       /* step 16 */  
   if (DEBUGLEVEL>2) fprintferr("Step 16\n");    if (DEBUGLEVEL>2) fprintferr("Step 16\n");
   K = FpM_ker(M, gell);    K = FpM_ker(M, gell);
   dK= lg(K)-1; if (!dK) { avma=av; return gzero; }    /* step 18 & ff */
       /* step 17 */  
   if (DEBUGLEVEL>2) fprintferr("Step 17\n");  
   listmodsup=cgetg(lSml2+1,t_VEC);  
   listbidsup=cgetg(lSml2+1,t_VEC);  
   listellranksup=cgetg(lSml2+1,t_VECSMALL);  
   for (i=1; i<=lSml2; i++)  
   {  
     listmodsup[i]=(long)idealmul(nfz,(GEN)listprSp[i],(GEN)listmod[i]);  
     listbidsup[i]=(long)zidealstarinitgen(nfz,(GEN)listmodsup[i]);  
     listellranksup[i] = ellrank(gmael3(listbidsup,i,2,2), ell);  
   }  
   vecMsup=cgetg(lSml2+1,t_VEC);  
   for (i=1; i<=lSml2; i++)  
   {  
     Msup=cgetg(nbcol+1,t_MAT); vecMsup[i]=(long)Msup;  
     for (j=1; j<=dv; j++) Msup[j]=lmul(mginv,ideallogauxsup(i,(GEN)vecw[j]));  
     for (   ; j<=nbcol; j++) Msup[j]=(long)ideallogauxsup(i,(GEN)vecalphap[j-dv]);  
   }  
       /* step 18 */  
   if (DEBUGLEVEL>2) fprintferr("Step 18\n");    if (DEBUGLEVEL>2) fprintferr("Step 18\n");
   do    dK = lg(K)-1;
     y = cgetg(dK+1,t_VECSMALL);
     res = cgetg(1, t_VEC); /* in case all = 1 */
     while (dK)
   {    {
     y = cgetg(dK,t_VECSMALL);  
     for (i=1; i<dK; i++) y[i] = 0;      for (i=1; i<dK; i++) y[i] = 0;
         /* step 19 */      y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
     for(;;)      do
     {      { /* cf. algo 5.3.18 */
       X = (GEN)K[dK];        GEN be, P, X = FpV_red(gmul_mati_smallvec(K, y), gell);
       for (j=1; j<dK; j++) X = gadd(X, gmulsg(y[j],(GEN)K[j]));        if (ok_congruence(X, gell, lW, vecMsup))
       X = FpV_red(X, gell);  
       finalresult = testx(bnf,X,module,subgroup,vecMsup);  
       if (finalresult) return gerepilecopy(av, finalresult);  
         /* step 20,21,22 */  
       i = dK;  
       do  
       {        {
         i--; if (!i) goto DECREASE;          be = compute_beta(X, vecWB, gell, bnfz);
         if (i < dK-1) y[i+1] = 0;          P = compute_polrel(nfz, &T, be, g, ell);
         y[i]++;          if (DEBUGLEVEL>1) fprintferr("polrel(beta) = %Z\n", P);
       } while (y[i] >= ell);          if (!all && gegal(subgroup, rnfnormgroup(bnr, P))) return P; /* DONE */
     }          res = concatsp(res, P);
 DECREASE:        }
     dK--;      } while (increment(y, dK, ell));
       y[dK--] = 0;
   }    }
   while (dK);    if (all) return res;
   avma = av; return gzero;    return gzero; /* FAIL */
   }
   
   GEN
   rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
   {
     gpmem_t av = avma;
     return gerepilecopy(av, _rnfkummer(bnr, subgroup, all, prec));
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>