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

Annotation of OpenXM_contrib/pari/src/modules/kummer.c, Revision 1.1

1.1     ! maekawa     1: /*******************************************************************/
        !             2: /*                                                                 */
        !             3: /*                      KUMMER EXTENSIONS                          */
        !             4: /*                                                                 */
        !             5: /*******************************************************************/
        !             6: /* $Id: kummer.c,v 1.1.1.1 1999/09/16 13:48:13 karim Exp $ */
        !             7: #include "pari.h"
        !             8:
        !             9: static long rc,gc,ell,degK,degKz,m,d,ru,rv,lraycyc,vnf,dv,nbcol,lSml2;
        !            10: static long lSl,lSp,dK,lSl2,dc;
        !            11: static GEN module,matexpoteta1,bnf,nf,raycyc,polnf,phiell,powtaubet,vecMsup;
        !            12: static GEN unmodell,bnfz,nfz,U,uu,gell,cyc,gencyc,vecalpha,R,A1,A2,g,vselmer;
        !            13: static GEN bnrz,polrel,steinitzZk,listmod,listprSp,listbid,listunif,listellrank;
        !            14: static GEN listmodsup,listbidsup,listellranksup,vecw,Sm,Sml1,Sml2,Sl;
        !            15: static GEN ESml2,Sp,vecbetap,vecalphap,matP,gg,mginv;
        !            16:
        !            17: /* row vector B x matrix T : c_j=prod_i (b_i^T_ij) */
        !            18: static GEN
        !            19: groupproduct(GEN B, GEN T)
        !            20: {
        !            21:   long lB,lT,i,j;
        !            22:   GEN c,p1;
        !            23:
        !            24:   lB=lg(B)-1;
        !            25:   lT=lg(T)-1;
        !            26:   c=cgetg(lT+1,t_VEC);
        !            27:   for (j=1; j<=lT; j++)
        !            28:   {
        !            29:     p1=gun;
        !            30:     for (i=1; i<=lB; i++) p1=gmul(p1,gpui((GEN)B[i],gcoeff(T,i,j),0));
        !            31:     c[j]=(long)p1;
        !            32:   }
        !            33:   return c;
        !            34: }
        !            35:
        !            36: static GEN
        !            37: grouppows(GEN B, long ex)
        !            38: {
        !            39:   long lB,j;
        !            40:   GEN c;
        !            41:
        !            42:   lB=lg(B)-1;
        !            43:   c=cgetg(lB+1,t_VEC);
        !            44:   for (j=1; j<=lB; j++) c[j]=(long)gpuigs((GEN)B[j],ex);
        !            45:   return c;
        !            46: }
        !            47:
        !            48: static GEN
        !            49: get_listx(GEN arch,GEN msign,GEN munit,GEN vecmunit2,long ell,long lSp,long lSml2,long nbvunit,long r1)
        !            50: {
        !            51:   GEN kermat,p2,p3,p4,y,yinit,listx,unmodell = gmodulss(1,ell);
        !            52:   long i,j,k,cmpt,lker;
        !            53:
        !            54:   kermat=lift(ker(munit)); lker=lg(kermat)-1;
        !            55:   if (!lker) return cgetg(1,t_VEC);
        !            56:   yinit=cgetg(lker,t_VEC); y=cgetg(lker,t_VEC);
        !            57:   for (i=1; i<lker; i++)
        !            58:   {
        !            59:     y[i]=lgeti(3); affii(gzero,(GEN)y[i]);
        !            60:     yinit[i]=zero;
        !            61:   }
        !            62:   listx=cgetg(1,t_VEC);
        !            63:   for(;;)
        !            64:   {
        !            65:     p2=cgetg(2,t_VEC);
        !            66:     p3=(GEN)kermat[lker];
        !            67:     for (j=1; j<lker; j++)
        !            68:       p3=lift(gmul(unmodell,gadd(p3,gmul((GEN)y[j],(GEN)kermat[j]))));
        !            69:     cmpt=0; for (j=1; j<=lSp; j++) if (gcmp0((GEN)p3[j+nbvunit])) cmpt++;
        !            70:     if (!cmpt)
        !            71:     {
        !            72:       for (k=1; k<=lSml2; k++)
        !            73:       {
        !            74:         p4=lift(gmul((GEN)vecmunit2[k],p3));
        !            75:         j=1; while (j<lg(p4) && !signe(p4[j])) j++;
        !            76:         if (j>=lg(p4)) break;
        !            77:       }
        !            78:       if (k>lSml2)
        !            79:       {
        !            80:         p4=lift(gmul(msign,p3));
        !            81:         j=1; while (j<=r1 && gegal((GEN)p4[j],(GEN)arch[j])) j++;
        !            82:         if (j>r1) { p2[1]=(long)p3; listx=concatsp(listx,p2); }
        !            83:       }
        !            84:     }
        !            85:     if (!lSp)
        !            86:     {
        !            87:       j=1; while (j<lker && !signe(y[j])) j++;
        !            88:       if (j<lker && gcmp1((GEN)y[j]))
        !            89:       {
        !            90:         p3=gsub(p3,(GEN)kermat[lker]);
        !            91:         for (k=1; k<=lSml2; k++)
        !            92:         {
        !            93:           p4=lift(gmul((GEN)vecmunit2[k],p3));
        !            94:           j=1; while (j<lg(p4) && !signe(p4[j])) j++;
        !            95:           if (j>=lg(p4)) break;
        !            96:         }
        !            97:         if (k>lSml2)
        !            98:         {
        !            99:           p4=lift(gmul(msign,p3));
        !           100:           j=1; while (j<=r1 && gegal((GEN)p4[j],(GEN)arch[j])) j++;
        !           101:           if (j>r1) { p2[1]=(long)p3; listx=concatsp(listx,p2); }
        !           102:         }
        !           103:       }
        !           104:     }
        !           105:     i=lker;
        !           106:     do
        !           107:     {
        !           108:       i--; if (!i) return listx;
        !           109:
        !           110:       if (i<lker-1) affii((GEN)yinit[i+1],(GEN)y[i+1]);
        !           111:       gaddgsz((GEN)y[i],1,(GEN)y[i]);
        !           112:     }
        !           113:     while (cmpis((GEN)y[i],ell)>=0);
        !           114:   }
        !           115: }
        !           116:
        !           117: static GEN
        !           118: reducealpha(GEN nf,GEN x,GEN gell)
        !           119: /* etant donne un nombre algebrique x du corps nf -- non necessairement
        !           120:    entier -- calcule un entier algebrique de la forme x*g^gell */
        !           121: {
        !           122:   long tx=typ(x),i;
        !           123:   GEN den,fa,fac,ep,p1,y;
        !           124:
        !           125:   nf=checknf(nf);
        !           126:   if ((tx==t_POL)||(tx==t_POLMOD)) y=algtobasis(nf,x);
        !           127:   else { y=x; x=basistoalg(nf,y); }
        !           128:   den=denom(y);
        !           129:   if (gcmp1(den)) return x;
        !           130:   fa=decomp(den);fac=(GEN)fa[1];ep=(GEN)fa[2];
        !           131:   p1=gun;
        !           132:   for (i=1; i<lg(fac); i++)
        !           133:     p1=mulii(p1,gpui((GEN)fac[i],gceil(gdiv((GEN)ep[i],gell)),0));
        !           134:   return gmul(gpui(p1,gell,0),x);
        !           135: }
        !           136:
        !           137: /* si all!=0, donne toutes les equations correspondant a un sousgroupe
        !           138:    de determinant all (i.e. de degre all) */
        !           139: static GEN
        !           140: rnfkummersimple(GEN bnr, GEN subgroup, long all, long prec)
        !           141: {
        !           142:   long r1,degnf,ell,j,i;
        !           143:   long nbgenclK,lprm,lprl,lSml2,lSl,lSp,lastindex,nbvunit,nbvunit0;
        !           144:   long kmax,nbgenbid,lastbid,nbcol,nblig,k,llistx,e,vp,vd,v1;
        !           145:
        !           146:   GEN bnf,nf,p1,classgroup,cyclic,bid,module,ideal,arch;
        !           147:   GEN gell,p2,p3,p22;
        !           148:   GEN clK,cyclicK,genK,listgamma,listalpha,fa,prm,expom,Sm,Sml1,Sml2,Sl,ESml2;
        !           149:   GEN pp,p,prl,Sp,vecbeta,matexpo,unmodell,vunit,id,pr,vsml2,listalpha0;
        !           150:   GEN cycpro,cl,cycbid,munit,vecmunit2,msign,archartif,listx,listal,listgamma0;
        !           151:   GEN vecbeta0;
        !           152:
        !           153:   checkbnrgen(bnr); bnf=(GEN)bnr[1];
        !           154:   nf=(GEN)bnf[7]; r1=itos(gmael(nf,2,1));
        !           155:   polnf=(GEN)nf[1]; vnf=varn(polnf);
        !           156:   if (vnf==0) err(talker,"main variable in kummer must not be x");
        !           157:   p1=conductor(bnr,all ? gzero : subgroup,2,prec);
        !           158:   bnr=(GEN)p1[2]; if(!all) subgroup=(GEN)p1[3];
        !           159:   classgroup=(GEN)bnr[5];
        !           160:   cyclic=(GEN)classgroup[2];
        !           161:   bid=(GEN)bnr[2]; module=(GEN)bid[1];
        !           162:   ideal=(GEN)module[1]; arch=(GEN)module[2];
        !           163:   if (!all && gcmp0(subgroup)) subgroup = diagonal(cyclic);
        !           164:   degnf=lgef(polnf)-3;
        !           165:   gell = all? stoi(all): det(subgroup);
        !           166:   if (!isprime(gell))
        !           167:     err(impl,"kummer for non prime relative degree");
        !           168:   ell=itos(gell);
        !           169:   unmodell=gmodulss(1,ell);
        !           170:   clK=gmael(bnf,8,1); cyclicK=(GEN)clK[2]; genK=(GEN)clK[3];
        !           171:   nbgenclK=lg(genK)-1;
        !           172:   i=1; while (i<=nbgenclK && gdivise((GEN)cyclicK[i],gell)) i++;
        !           173:   lastindex=i-1;
        !           174:   listgamma0=cgetg(nbgenclK+1,t_VEC);
        !           175:   listgamma=cgetg(nbgenclK+1,t_VEC);
        !           176:   listalpha0=cgetg(lastindex+1,t_VEC);
        !           177:   listalpha=cgetg(lastindex+1,t_VEC);
        !           178:   p1=gmul(gell,ideal);
        !           179:   for (i=1; i<=lastindex; i++)
        !           180:   {
        !           181:     p3=basistoalg(nf,idealcoprime(nf,(GEN)genK[i],p1));
        !           182:     listgamma[i]=listgamma0[i]=linv(p3);
        !           183:     p2=idealpow(nf,(GEN)genK[i],(GEN)cyclicK[i]);
        !           184:     listalpha0[i]=(long)basistoalg(nf,(GEN)isprincipalgen(bnf,p2)[2]);
        !           185:     listalpha[i]=lmul((GEN)listalpha0[i],gpui(p3,(GEN)cyclicK[i],0));
        !           186:   }
        !           187:   for (   ; i<=nbgenclK; i++)
        !           188:   {
        !           189:     p3=basistoalg(nf,idealcoprime(nf,(GEN)genK[i],p1));
        !           190:     p2=idealpow(nf,(GEN)genK[i],(GEN)cyclicK[i]);
        !           191:     p2=basistoalg(nf,(GEN)isprincipalgen(bnf,p2)[2]);
        !           192:     p22=lift(ginv(gmul(unmodell,(GEN)cyclicK[i])));
        !           193:     p2=gpui(p2,p22,0);listgamma0[i]=(long)p2;
        !           194:     listgamma[i]=lmul(p2,gpui(p3,addis(mulii(p22,(GEN)cyclicK[i]),-1),0));
        !           195:   }
        !           196:   /* a ce stade, module est egal au conducteur. */
        !           197:   fa=idealfactor(nf,ideal);
        !           198:   prm=(GEN)fa[1]; expom=(GEN)fa[2]; lprm=lg(prm)-1;
        !           199:   Sm=cgetg(1,t_VEC);Sml1=cgetg(1,t_VEC);Sml2=cgetg(1,t_VEC);Sl=cgetg(1,t_VEC);
        !           200:   ESml2=cgetg(1,t_VEC);
        !           201:   for (i=1; i<=lprm; i++)
        !           202:   {
        !           203:     pp=cgetg(2,t_VEC); pp[1]=prm[i];
        !           204:     p=gmael(pp,1,1); e=itos(gmael(pp,1,3));
        !           205:     if (!gegal(p,gell))
        !           206:     {
        !           207:       if (!gcmp1((GEN)expom[i]))
        !           208:       {
        !           209:        if (all) return cgetg(1,t_VEC);
        !           210:        err(talker,"bug1 in kummer");
        !           211:       }
        !           212:       Sm=concatsp(Sm,pp);
        !           213:     }
        !           214:     else
        !           215:     {
        !           216:       vp=itos((GEN)expom[i]); vd=(vp-1)*(ell-1)-ell*e;
        !           217:       if (vd > 0)
        !           218:       {
        !           219:        if (all) return cgetg(1,t_VEC);
        !           220:        err(talker,"bug4 in kummer");
        !           221:       }
        !           222:       if (vd==0)
        !           223:         Sml1=concatsp(Sml1,pp);
        !           224:       else
        !           225:       {
        !           226:        if (vp==1)
        !           227:        {
        !           228:          if (all) return cgetg(1,t_VEC);
        !           229:          err(talker,"bug2 in kummer");
        !           230:        }
        !           231:        else
        !           232:        {
        !           233:          Sml2=concatsp(Sml2,pp); ESml2=concatsp(ESml2,(GEN)expom[i]);
        !           234:        }
        !           235:       }
        !           236:     }
        !           237:   }
        !           238:   prl=primedec(nf,gell); lprl=lg(prl)-1;
        !           239:   for (i=1; i<=lprl; i++)
        !           240:     if (!idealval(nf,ideal,(GEN)prl[i]))
        !           241:     {
        !           242:       pp=cgetg(2,t_VEC); pp[1]=prl[i];
        !           243:       Sl=concatsp(Sl,pp);
        !           244:     }
        !           245:   lSml2=lg(Sml2)-1; lSl=lg(Sl)-1;
        !           246:   Sp=concatsp(Sm,Sml1); lSp=lg(Sp)-1;
        !           247:   vecbeta=cgetg(lSp+1,t_VEC); matexpo=cgetg(lSp+1,t_MAT);
        !           248:   vecbeta0=cgetg(lSp+1,t_VEC);
        !           249:   for (j=1; j<=lSp; j++)
        !           250:   {
        !           251:     p1=isprincipalgen(bnf,(GEN)Sp[j]);
        !           252:     p2=basistoalg(nf,(GEN)p1[2]);
        !           253:     for (i=1; i<=lastindex; i++)
        !           254:       p2 = gmul(p2,gpui((GEN)listgamma[i],gmael(p1,1,i),0));
        !           255:     p22=p2;
        !           256:     for (   ; i<=nbgenclK; i++)
        !           257:     {
        !           258:       p2 = gmul(p2,gpui((GEN)listgamma[i],gmael(p1,1,i),0));
        !           259:       p22 = gmul(p22,gpui((GEN)listgamma0[i],gmael(p1,1,i),0));
        !           260:     }
        !           261:     matexpo[j]=(long)p1[1];
        !           262:     vecbeta[j]=(long)p2; /* attention, ceci sont les beta modifies */
        !           263:     vecbeta0[j]=(long)p22;
        !           264:   }
        !           265:   matexpo=gmul(unmodell,matexpo);
        !           266:   vunit=gmodulcp(concatsp(gmael(bnf,8,5),gmael3(bnf,8,4,2)),polnf);
        !           267:   nbvunit0=lg(vunit)-1;
        !           268:   for (i=1; i<=lastindex; i++) vunit=concatsp(vunit,(GEN)listalpha[i]);
        !           269:   nbvunit=lg(vunit)-1;
        !           270:   id=idmat(degnf);
        !           271:   for (i=1; i<=lSl; i++)
        !           272:   {
        !           273:     pr=(GEN)Sl[i];
        !           274:     id=idealmul(nf,idealpow(nf,pr,stoi((ell*itos((GEN)pr[3]))/(ell-1))),id);
        !           275:   }
        !           276:   vsml2=cgetg(lSml2+1,t_VEC);
        !           277:   for (i=1; i<=lSml2; i++)
        !           278:   {
        !           279:     pr=(GEN)Sml2[i];
        !           280:     v1=itos((GEN)ESml2[i]);
        !           281:     kmax=(itos((GEN)pr[3])*ell)/(ell-1)+1-v1;
        !           282:     p1=idealpow(nf,pr,stoi(kmax));
        !           283:     id=idealmul(nf,id,p1);
        !           284:     p2=idealmul(nf,p1,pr);
        !           285:     p3=zidealstarinitgen(nf,p2);
        !           286:     vsml2[i]=(long)p3;
        !           287:     cycpro=gmael(p3,2,2);
        !           288:     for (j=1;j<lg(cycpro);j++)
        !           289:       if (! gdivise((GEN)cycpro[j],gell)) err(talker,"bug5 in kummer");
        !           290:   }
        !           291:   bid=zidealstarinitgen(nf,id);
        !           292:   cl=(GEN)bid[2]; cycbid=(GEN)cl[2]; nbgenbid=lg(cycbid)-1;
        !           293:   i=1; while (i<=nbgenbid && gdivise((GEN)cycbid[i],gell)) i++;
        !           294:   lastbid=i-1;
        !           295:   nbcol=nbvunit+lSp; nblig=lastbid+lastindex;
        !           296:   munit=cgetg(nbcol+1,t_MAT); vecmunit2=cgetg(lSml2+1,t_VEC);
        !           297:   msign=cgetg(nbcol+1,t_MAT);
        !           298:   for (k=1; k<=lSml2; k++) vecmunit2[k]=lgetg(nbcol+1,t_MAT);
        !           299:   archartif=cgetg(r1+1,t_VEC); for (j=1; j<=r1; j++) archartif[j]=un;
        !           300:   for (j=1; j<=nbvunit; j++)
        !           301:   {
        !           302:     p1=zideallog(nf,(GEN)vunit[j],bid);
        !           303:     p2=cgetg(nblig+1,t_COL); munit[j]=(long)p2;
        !           304:     for (i=1; i<=lastbid; i++) p2[i]=p1[i];
        !           305:     for (   ; i<=nblig; i++) p2[i]=zero;
        !           306:     for (k=1; k<=lSml2; k++)
        !           307:       ((GEN)vecmunit2[k])[j]=(long)zideallog(nf,(GEN)vunit[j],(GEN)vsml2[k]);
        !           308:     msign[j]=(long)zsigne(nf,(GEN)vunit[j],archartif);
        !           309:   }
        !           310:   p1=cgetg(nbgenclK-lastindex+1,t_VEC);
        !           311:   for (j=1; j<=lSp; j++)
        !           312:   {
        !           313:     p2=zideallog(nf,(GEN)vecbeta[j],bid);
        !           314:     p3=cgetg(nblig+1,t_COL); munit[j+nbvunit]=(long)p3;
        !           315:     for (i=1; i<=lastbid; i++) p3[i]=p2[i];
        !           316:     for (   ; i<=nblig; i++) p3[i]=coeff(matexpo,i-lastbid,j);
        !           317:     for (k=1; k<=lSml2; k++)
        !           318:     {
        !           319:       p2=zideallog(nf,(GEN)vecbeta[j],(GEN)vsml2[k]);
        !           320:       ((GEN)vecmunit2[k])[j+nbvunit]=(long)p2;
        !           321:     }
        !           322:     msign[j+nbvunit]=(long)zsigne(nf,(GEN)vecbeta[j],archartif);
        !           323:   }
        !           324:   munit=gmul(unmodell,munit); vecmunit2=gmul(unmodell,vecmunit2);
        !           325:   listx = get_listx(arch,msign,munit,vecmunit2,ell,lSp,lSml2,nbvunit,r1);
        !           326:   llistx=lg(listx)-1;
        !           327:   listal=cgetg(llistx+1,t_VEC);
        !           328:   for (i=1; i<=llistx; i++)
        !           329:   {
        !           330:     p1=gun; p3=(GEN)listx[i];
        !           331:     for (j=1; j<=nbvunit0; j++)
        !           332:       p1 = gmul(p1,gpui((GEN)vunit[j],(GEN)p3[j],0));
        !           333:     for (   ; j<=nbvunit; j++)
        !           334:       p1 = gmul(p1,gpui((GEN)listalpha0[j-nbvunit0],(GEN)p3[j],0));
        !           335:     for (j=1; j<=lSp; j++)
        !           336:       p1 = gmul(p1,gpui((GEN)vecbeta0[j],(GEN)p3[j+nbvunit],0));
        !           337:     listal[i]=(long)reducealpha(nf,p1,gell);
        !           338:   }
        !           339:  /* A ce stade, tous les alpha dans la liste listal statisfont les
        !           340:   * congruences, noncongruences et signatures, et ne sont pas des puissances
        !           341:   * l-iemes donc x^l-alpha est irreductible, sa signature est correcte ainsi
        !           342:   * que son discriminant relatif. Reste a determiner son groupe de normes.
        !           343:   * Les alpha ne sont peut-etre pas des entiers algebriques..  */
        !           344:   if (DEBUGLEVEL)
        !           345:     { fprintferr("listalpha = \n"); outerr(listal); flusherr(); }
        !           346:   p2=cgetg(1,t_VEC);
        !           347:   for (i=1; i<=llistx; i++)
        !           348:   {
        !           349:     p1=gsub(gpuigs(polx[0],ell),(GEN)listal[i]);
        !           350:     if (all || gegal(rnfnormgroup(bnr,p1),subgroup))
        !           351:       p2=concatsp(p2,p1);
        !           352:   }
        !           353:   if (all) return gcopy(p2);
        !           354:   switch(lg(p2)-1)
        !           355:   {
        !           356:     case 0: err(talker,"bug 6: no equation found in kummer");
        !           357:     case 1: return gcopy((GEN)p2[1]);
        !           358:     default: fprintferr("equations = \n"); outerr(p2);
        !           359:       err(talker,"bug 7: more than one equation found in kummer");
        !           360:   }
        !           361:   return gzero; /* not reached */
        !           362: }
        !           363:
        !           364: static GEN
        !           365: tauofideal(GEN id)
        !           366: {
        !           367:   long j;
        !           368:   GEN p1,p2;
        !           369:
        !           370:   p1=gsubst(gmul((GEN)nfz[7],id),vnf,U);
        !           371:   p2=cgetg(lg(p1),t_MAT);
        !           372:   for (j=1; j<lg(p1); j++) p2[j]=(long)algtobasis(nfz,(GEN)p1[j]);
        !           373:   return p2;
        !           374: }
        !           375:
        !           376: static GEN
        !           377: tauofprimeideal(GEN pr)
        !           378: {
        !           379:   GEN p1;
        !           380:
        !           381:   p1=gcopy(pr);
        !           382:   p1[2]=(long)algtobasis(nfz,gsubst(gmul((GEN)nfz[7],(GEN)pr[2]),vnf,U));
        !           383:   return gcoeff(idealfactor(nfz,prime_to_ideal(nfz,p1)),1,1);
        !           384: }
        !           385:
        !           386: static long
        !           387: isprimeidealconj(GEN pr1, GEN pr2)
        !           388: {
        !           389:   GEN pr=pr1;
        !           390:
        !           391:   do
        !           392:   {
        !           393:     if (gegal(pr,pr2)) return 1;
        !           394:     pr=tauofprimeideal(pr);
        !           395:   }
        !           396:   while (!gegal(pr,pr1));
        !           397:   return 0;
        !           398: }
        !           399:
        !           400: static long
        !           401: isconjinprimelist(GEN listpr, GEN pr2)
        !           402: {
        !           403:   long ll=lg(listpr)-1,i;
        !           404:
        !           405:   for (i=1; i<=ll; i++)
        !           406:     if (isprimeidealconj((GEN)listpr[i],pr2)) return 1;
        !           407:   return 0;
        !           408: }
        !           409:
        !           410: static void
        !           411: computematexpoteta1(void)
        !           412: {
        !           413:   long i,j;
        !           414:   GEN p4,p5,p6;
        !           415:
        !           416:   matexpoteta1=cgetg(degK+1,t_MAT);
        !           417:   p4=gmodulcp(polun[vnf],R);
        !           418:   for (j=1; j<=degK; j++)
        !           419:   {
        !           420:     p5=cgetg(degKz+1,t_COL); matexpoteta1[j]=(long)p5; p6=lift(p4);
        !           421:     for (i=1; i<=degKz; i++) p5[i]=(long)truecoeff(p6,i-1);
        !           422:     if (j<degK) p4=gmul(p4,A1);
        !           423:   }
        !           424: }
        !           425:
        !           426: static GEN
        !           427: downtoK(GEN x)
        !           428: {
        !           429:   long i;
        !           430:   GEN p1,p2,p3;
        !           431:
        !           432:   p1=cgetg(degKz+1,t_COL);
        !           433:   p2=lift(x);
        !           434:   for (i=1; i<=degKz; i++) p1[i]=(long)truecoeff(p2,i-1);
        !           435:   p2=inverseimage(matexpoteta1,p1);
        !           436:   if (lg(p2)==1) err(talker,"not an element of K in downtoK");
        !           437:   p3=(GEN)p2[degK];
        !           438:   for (i=degK-1; i; i--) p3=gadd((GEN)p2[i],gmul(polx[vnf],p3));
        !           439:   return gmodulcp(p3,polnf);
        !           440: }
        !           441:
        !           442: static GEN
        !           443: tracetoK(GEN x)
        !           444: {
        !           445:   long i;
        !           446:   GEN p1,p2;
        !           447:
        !           448:   p1=x; p2=x;
        !           449:   for (i=1; i<=m-1; i++)
        !           450:   {
        !           451:     p1=gsubst(lift(p1),vnf,U);
        !           452:     p2=gadd(p2,p1);
        !           453:   }
        !           454:   return downtoK(p2);
        !           455: }
        !           456:
        !           457: static GEN
        !           458: normtoK(GEN x)
        !           459: {
        !           460:   long i;
        !           461:   GEN p1,p2;
        !           462:
        !           463:   p1=x; p2=x;
        !           464:   for (i=1; i<=m-1; i++)
        !           465:   {
        !           466:     p1=gsubst(lift(p1),vnf,U);
        !           467:     p2=gmul(p2,p1);
        !           468:   }
        !           469:   return downtoK(p2);
        !           470: }
        !           471:
        !           472: static GEN
        !           473: computepolrel(void)
        !           474: {
        !           475:   long i,j;
        !           476:   GEN p1,p2;
        !           477:
        !           478:   p1=gun; p2=gmodulcp(polx[vnf],R);
        !           479:   for (i=0; i<=m-1; i++)
        !           480:   {
        !           481:     p1=gmul(p1,gsub(polx[0],p2));
        !           482:     if (i<m-1) p2=gsubst(lift(p2),vnf,U);
        !           483:   }
        !           484:   for (j=2; j<=m+2; j++) p1[j]=(long)downtoK((GEN)p1[j]);
        !           485:   return p1;
        !           486: }
        !           487:
        !           488: /* alg. 5.2.15. with remark */
        !           489: static GEN
        !           490: isprincipalell(GEN id)
        !           491: {
        !           492:   long i;
        !           493:   GEN p1,p2,logdisc,be;
        !           494:
        !           495:   p1=isprincipalgen(bnfz,id); logdisc=(GEN)p1[1];
        !           496:   be=basistoalg(bnfz,(GEN)p1[2]);
        !           497:   for (i=rc+1; i<=gc; i++)
        !           498:     be=gmul(be,gpui((GEN)vecalpha[i],modii(mulii((GEN)logdisc[i],(GEN)uu[i]),gell),0));
        !           499:   p2=cgetg(3,t_VEC);
        !           500:   p2[2]=(long)be;
        !           501:   p1=cgetg(rc+1,t_COL); p2[1]=(long)p1;
        !           502:   for (i=1; i<=rc; i++) p1[i]=logdisc[i];
        !           503:   return p2;
        !           504: }
        !           505:
        !           506: /* alg. 5.3.11. */
        !           507: static GEN
        !           508: isvirtualunit(GEN v)
        !           509: {
        !           510:   long llist,i,j,ex;
        !           511:   GEN p1,p2,listex,listpr,q,ga,eps,vecy,logdisc;
        !           512:
        !           513:   p1=idealfactor(nfz,v);
        !           514:   listex=(GEN)p1[2]; listpr=(GEN)p1[1]; llist=lg(listex)-1;
        !           515:   q=idmat(degKz);
        !           516:   for (i=1; i<=llist; i++)
        !           517:   {
        !           518:     ex=itos((GEN)listex[i]);
        !           519:     if (ex%ell) err(talker,"not a virtual unit in isvirtualunit");
        !           520:     q=idealmul(nfz,q,idealpow(nfz,(GEN)listpr[i],stoi(ex/ell)));
        !           521:   }
        !           522:   p1=isprincipalgen(bnfz,q); logdisc=(GEN)p1[1];
        !           523:   ga=basistoalg(nfz,(GEN)p1[2]);
        !           524:   for (j=rc+1; j<=gc; j++)
        !           525:     ga=gmul(ga,gpui((GEN)vecalpha[j],divii((GEN)logdisc[j],(GEN)cyc[j]),0));
        !           526:   eps=gpuigs(ga,ell);
        !           527:   vecy=cgetg(rc+1,t_COL);
        !           528:   for (j=1; j<=rc; j++)
        !           529:   {
        !           530:     vecy[j]=(long)divii((GEN)logdisc[j],divii((GEN)cyc[j],gell));
        !           531:     eps=gmul(eps,gpui((GEN)vecalpha[j],(GEN)vecy[j],0));
        !           532:   }
        !           533:   eps=gdiv(v,eps);
        !           534:   p1=isunit(bnfz,eps);
        !           535:   p2=cgetg(3,t_VEC);
        !           536:   p2[1]=(long)concatsp(vecy,lift(p1));
        !           537:   p2[2]=(long)ga;
        !           538:   return p2;
        !           539: }
        !           540:
        !           541: static GEN
        !           542: lifttokz(GEN id)
        !           543: {
        !           544:   long i,j;
        !           545:   GEN p1,p2,p3;
        !           546:
        !           547:   p1=gsubst(gmul((GEN)nf[7],id),vnf,A1);
        !           548:   p2=gmodulcp((GEN)nfz[7],R);
        !           549:   p3=cgetg(degK*degKz+1,t_MAT);
        !           550:   for (i=1; i<=degK; i++)
        !           551:     for (j=1; j<=degKz; j++)
        !           552:       p3[(i-1)*degKz+j]=(long)algtobasis(nfz,gmul((GEN)p1[i],(GEN)p2[j]));
        !           553:   return hnfmod(p3,detint(p3));
        !           554: }
        !           555:
        !           556: static GEN
        !           557: steinitzaux(GEN id)
        !           558: {
        !           559:   long i,j;
        !           560:   GEN p1,p2,p3,vecid,matid,pseudomat,pid;
        !           561:
        !           562:   p1=gsubst(gmul((GEN)nfz[7],id),vnf,polx[0]);
        !           563:   for (j=1; j<=degKz; j++)
        !           564:     p1[j]=(long)gmod((GEN)p1[j],polrel);
        !           565:   p2=cgetg(degKz+1,t_MAT);
        !           566:   for (j=1; j<=degKz; j++)
        !           567:   {
        !           568:     p3=cgetg(m+1,t_COL); p2[j]=(long)p3;
        !           569:     for (i=1; i<=m; i++) p3[i]=(long)algtobasis(nf,truecoeff((GEN)p1[j],i-1));
        !           570:   }
        !           571:   vecid=cgetg(degKz+1,t_VEC); matid=idmat(degK);
        !           572:   for (j=1; j<=degKz; j++) vecid[j]=(long)matid;
        !           573:   pseudomat=cgetg(3,t_VEC);
        !           574:   pseudomat[1]=(long)p2; pseudomat[2]=(long)vecid;
        !           575:   pid=(GEN)nfhermite(nf,pseudomat)[2];
        !           576:   p1=matid;
        !           577:   for (j=1; j<=m; j++) p1=idealmul(nf,p1,(GEN)pid[j]);
        !           578:   return p1;
        !           579: }
        !           580:
        !           581: static GEN
        !           582: normrelz(GEN id)
        !           583: {
        !           584:   GEN p1;
        !           585:
        !           586:   p1=steinitzaux(id);
        !           587:   return idealdiv(nf,p1,steinitzZk);
        !           588: }
        !           589:
        !           590: static GEN
        !           591: invimsubgroup(GEN bnr, GEN subgroup, GEN idealz, long prec)
        !           592: {
        !           593:   long lraycycz,i,j;
        !           594:   GEN Plog,rayclgpz,genraycycz,utemp,p1,p2;
        !           595:
        !           596:   bnrz=buchrayinitgen(bnfz,idealz,prec);
        !           597:   rayclgpz=(GEN)bnrz[5];
        !           598:   genraycycz=(GEN)rayclgpz[3]; lraycycz=lg(genraycycz)-1;
        !           599:   Plog=cgetg(lraycycz+lraycyc+1,t_MAT);
        !           600:   for (j=1; j<=lraycycz; j++)
        !           601:     Plog[j]=(long)isprincipalray(bnr,normrelz((GEN)genraycycz[j]));
        !           602:   for (   ; j<=lraycycz+lraycyc; j++)
        !           603:     Plog[j]=subgroup[j-lraycycz];
        !           604:   utemp=(GEN)hnfall(Plog)[2];
        !           605:   p2=cgetg((lraycycz<<1)+1,t_MAT);
        !           606:   for (j=1; j<=lraycycz; j++)
        !           607:   {
        !           608:     p1=cgetg(lraycycz+1,t_COL); p2[j]=(long)p1;
        !           609:     for (i=1; i<=lraycycz; i++) p1[i]=coeff(utemp,i,j);
        !           610:   }
        !           611:   p1=diagonal((GEN)rayclgpz[2]);
        !           612:   for (  ; j<=(lraycycz<<1); j++) p2[j]=p1[j-lraycycz];
        !           613:   return hnfmod(p2,(GEN)rayclgpz[1]);
        !           614: }
        !           615:
        !           616: static GEN
        !           617: ideallogaux(long ind, GEN al)
        !           618: {
        !           619:   long llogli,valal;
        !           620:   GEN p1;
        !           621:
        !           622:   valal=element_val(nfz,algtobasis(nfz,al),(GEN)listprSp[ind]);
        !           623:   al=gmul(al,gpuigs((GEN)listunif[ind],valal));
        !           624:   p1=zideallog(nfz,al,(GEN)listbid[ind]);
        !           625:   llogli=itos((GEN)listellrank[ind]);
        !           626:   setlg(p1,llogli+1); return p1;
        !           627: }
        !           628:
        !           629: static GEN
        !           630: ideallogauxsup(long ind, GEN al)
        !           631: {
        !           632:   long llogli,valal;
        !           633:   GEN p1;
        !           634:
        !           635:   valal=element_val(nfz,algtobasis(nfz,al),(GEN)listprSp[ind]);
        !           636:   al=gmul(al,gpuigs((GEN)listunif[ind],valal));
        !           637:   p1=zideallog(nfz,al,(GEN)listbidsup[ind]);
        !           638:   llogli=itos((GEN)listellranksup[ind]);
        !           639:   setlg(p1,llogli+1); return p1;
        !           640: }
        !           641:
        !           642: static GEN
        !           643: vectau(GEN list)
        !           644: {
        !           645:   long i,j,n,cmpt;
        !           646:   GEN listz,listc,yz,yc,listfl,s,p1,p2,y;
        !           647:
        !           648:   listz=(GEN)list[1]; listc=(GEN)list[2]; n=lg(listz)-1;
        !           649:   yz=cgetg(n+1,t_VEC); yc=cgetg(n+1,t_VEC);
        !           650:   listfl=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) listfl[i]=un;
        !           651:   cmpt=0;
        !           652:   for (i=1; i<=n; i++)
        !           653:   {
        !           654:     if (signe((GEN)listfl[i]))
        !           655:     {
        !           656:       cmpt++;
        !           657:       yz[cmpt]=listz[i];
        !           658:       s=(GEN)listc[i];
        !           659:       for (j=i+1; j<=n; j++)
        !           660:       {
        !           661:        if (signe((GEN)listfl[j])&&gegal((GEN)listz[j],(GEN)listz[i]))
        !           662:        {
        !           663:          s=gadd(s,(GEN)listc[j]);
        !           664:          listfl[j]=zero;
        !           665:        }
        !           666:       }
        !           667:       yc[cmpt]=(long)s;
        !           668:     }
        !           669:   }
        !           670:   y=cgetg(3,t_VEC);
        !           671:   p1=cgetg(cmpt+1,t_VEC); p2=cgetg(cmpt+1,t_VEC); y[1]=(long)p1; y[2]=(long)p2;
        !           672:   for (i=1; i<=cmpt; i++)
        !           673:   {
        !           674:     p1[i]=yz[i];
        !           675:     p2[i]=yc[i];
        !           676:   }
        !           677:   return y;
        !           678: }
        !           679:
        !           680: #if 0
        !           681: static GEN
        !           682: addtau(GEN listx, GEN listy)
        !           683: {
        !           684:   GEN p1,p2,y;
        !           685:
        !           686:   p1=concat((GEN)listx[1],(GEN)listy[1]);
        !           687:   p2=concat((GEN)listx[2],(GEN)listy[2]);
        !           688:   y=cgetg(3,t_VEC); y[1]=(long)p1; y[2]=(long)p2;
        !           689:   return vectau(y);
        !           690: }
        !           691: #endif
        !           692:
        !           693: static GEN
        !           694: subtau(GEN listx, GEN listy)
        !           695: {
        !           696:   GEN p1,p2,y;
        !           697:
        !           698:   p1=concat((GEN)listx[1],(GEN)listy[1]);
        !           699:   p2=concat((GEN)listx[2],gneg_i((GEN)listy[2]));
        !           700:   y=cgetg(3,t_VEC); y[1]=(long)p1; y[2]=(long)p2;
        !           701:   return vectau(y);
        !           702: }
        !           703:
        !           704: static GEN
        !           705: negtau(GEN list)
        !           706: {
        !           707:   GEN y;
        !           708:
        !           709:   y=cgetg(3,t_VEC);
        !           710:   y[1]=list[1];
        !           711:   y[2]=lneg((GEN)list[2]);
        !           712:   return y;
        !           713: }
        !           714:
        !           715: static GEN
        !           716: multau(GEN listx, GEN listy)
        !           717: {
        !           718:   GEN lzx,lzy,lcx,lcy,lzz,lcz,y;
        !           719:   long nx,ny,i,j,k;
        !           720:
        !           721:   lzx=(GEN)listx[1]; lcx=(GEN)listx[2]; nx=lg(lzx)-1;
        !           722:   lzy=(GEN)listy[1]; lcy=(GEN)listy[2]; ny=lg(lzy)-1;
        !           723:   lzz=cgetg(nx*ny+1,t_VEC); lcz=cgetg(nx*ny+1,t_VEC);
        !           724:   for (i=1; i<=nx; i++)
        !           725:     for (j=1; j<=ny; j++)
        !           726:     {
        !           727:       k=(i-1)*ny+j;
        !           728:       lzz[k]=ladd((GEN)lzx[i],(GEN)lzy[j]);
        !           729:       lcz[k]=lmul((GEN)lcx[i],(GEN)lcy[j]);
        !           730:     }
        !           731:   y=cgetg(3,t_VEC); y[1]=(long)lzz; y[2]=(long)lcz;
        !           732:   return vectau(y);
        !           733: }
        !           734:
        !           735: #if 0
        !           736: static GEN
        !           737: puipol(GEN be, GEN pol)
        !           738: {
        !           739:   long i,degpol;
        !           740:   GEN y;
        !           741:
        !           742:   degpol=degree(pol);
        !           743:   y=gun;
        !           744:   for (i=0; i<=degpol; i++)
        !           745:     if (!gcmp0((GEN)pol[i+2])) y=gmul(y,gpui((GEN)powtaubet[i+1],(GEN)pol[i+2],0));
        !           746:   return y;
        !           747: }
        !           748:
        !           749: static GEN
        !           750: puitau(GEN be, GEN list)
        !           751: {
        !           752:   long lgl,i;
        !           753:   GEN y,listz,listc,p1;
        !           754:
        !           755:   y=gun; listz=(GEN)list[1]; listc=(GEN)list[2]; lgl=lg(listz)-1;
        !           756:   for (i=1; i<=lgl; i++)
        !           757:   {
        !           758:     p1=lift((GEN)listc[i]);
        !           759:     if (!gcmp0(p1))
        !           760:     {
        !           761:       if (!gcmp0(lift(gmul(unmodell,(GEN)listz[i]))))
        !           762:        err(talker,"bug1 in puitau");
        !           763:       if (degree(p1)) err(talker,"bug2 in puitau");
        !           764:       if (typ(p1)==t_POL) p1=(GEN)p1[2];
        !           765:       y=gmul(y,gpui(puipol(be,(GEN)listz[i]),p1,0));
        !           766:     }
        !           767:   }
        !           768:   return y;
        !           769: }
        !           770: #endif
        !           771:
        !           772: static GEN
        !           773: mulpoltau(GEN poltau, GEN list)
        !           774: {
        !           775:   long i,j;
        !           776:   GEN y;
        !           777:
        !           778:   j=lg(poltau)-2;
        !           779:   y=cgetg(j+3,t_VEC);
        !           780:   y[1]=(long)negtau(multau(list,(GEN)poltau[1]));
        !           781:   for (i=2; i<=j+1; i++)
        !           782:   {
        !           783:     y[i]=(long)subtau((GEN)poltau[i-1],multau(list,(GEN)poltau[i]));
        !           784:   }
        !           785:   y[j+2]=poltau[j+1];
        !           786:   return y;
        !           787: }
        !           788:
        !           789: /* th. 5.3.5. and prop. 5.3.9. */
        !           790: static GEN
        !           791: computepolrelbeta(GEN be)
        !           792: {
        !           793:   long i,a,b,j;
        !           794:   GEN e,u,u1,u2,u3,p1,p2,p3,p4,zet,be1,be2,listr,s,veczi,vecci;
        !           795:
        !           796:   switch (ell)
        !           797:   {
        !           798:     case 2: err(talker,"you should not be here in rnfkummer !!"); break;
        !           799:     case 3: e=normtoK(be); u=tracetoK(be);
        !           800:       return gsub(gmul(polx[0],gsub(gsqr(polx[0]),gmulsg(3,e))),gmul(e,u));
        !           801:     case 5: e=normtoK(be);
        !           802:       if (d==2)
        !           803:       {
        !           804:        u=tracetoK(gpuigs(be,3));
        !           805:        p1=gadd(gmulsg(5,gsqr(e)),gmul(gsqr(polx[0]),gsub(gsqr(polx[0]),gmulsg(5,e))));
        !           806:        return gsub(gmul(polx[0],p1),gmul(e,u));
        !           807:       }
        !           808:       else
        !           809:       {
        !           810:        be1=gsubst(lift(be),vnf,U); be2=gsubst(lift(be1),vnf,U);
        !           811:        u1=tracetoK(gmul(be,be1)); u2=tracetoK(gmul(gmul(be,be2),gsqr(be1)));
        !           812:        u3=tracetoK(gmul(gmul(gsqr(be),gpuigs(be1,3)),be2));
        !           813:        p1=gsub(gsqr(polx[0]),gmulsg(10,e));
        !           814:        p1=gsub(gmul(polx[0],p1),gmulsg(5,gmul(e,u1)));
        !           815:        p1=gadd(gmul(polx[0],p1),gmul(gmulsg(5,e),gsub(e,u2)));
        !           816:        p1=gsub(gmul(polx[0],p1),gmul(e,u3));
        !           817:        return p1;
        !           818:       }
        !           819:     default: p1=cgetg(2,t_VEC); p2=cgetg(3,t_VEC); p3=cgetg(2,t_VEC);
        !           820:       p4=cgetg(2,t_VEC); p3[1]=zero; p4[1]=un;
        !           821:       p2[1]=(long)p3; p2[2]=(long)p4; p1[1]=(long)p2;
        !           822:       zet=gmodulcp(polx[vnf],phiell);
        !           823:       listr=cgetg(m+1,t_VEC);
        !           824:       listr[1]=un;
        !           825:       for (i=2; i<=m; i++) listr[i]=(long)modii(mulii((GEN)listr[i-1],g),gell);
        !           826:       veczi=cgetg(m+1,t_VEC);
        !           827:       for (b=0; b<m; b++)
        !           828:       {
        !           829:        s=gzero;
        !           830:        for (a=0; a<m; a++)
        !           831:          s=gadd(gmul(polx[0],s),modii(mulii((GEN)listr[b+1],(GEN)listr[a+1]),gell));
        !           832:        veczi[b+1]=(long)s;
        !           833:       }
        !           834:       for (j=0; j<ell; j++)
        !           835:       {
        !           836:        vecci=cgetg(m+1,t_VEC);
        !           837:        for (b=0; b<m; b++) vecci[b+1]=(long)gpui(zet,mulsi(j,(GEN)listr[b+1]),0);
        !           838:        p4=cgetg(3,t_VEC);
        !           839:        p4[1]=(long)veczi; p4[2]=(long)vecci;
        !           840:        p1=mulpoltau(p1,p4);
        !           841:       }
        !           842:       powtaubet=cgetg(m+1,t_VEC);
        !           843:       powtaubet[1]=(long)be;
        !           844:       for (i=2; i<=m; i++) powtaubet[i]=(long)gsubst(lift((GEN)powtaubet[i-1]),vnf,U);
        !           845:       err(impl,"difficult Kummer for ell>=7"); return gzero;
        !           846:   }
        !           847:   return NULL; /* not reached */
        !           848: }
        !           849:
        !           850: /* multiply be by ell-th powers of units as to find "small" L2-norm for new be */
        !           851: static GEN
        !           852: reducebeta(GEN be)
        !           853: {
        !           854:   long i,lu;
        !           855:   GEN p1,p2,unitsell,unitsellinv,nmax,ben;
        !           856:
        !           857:   unitsell=grouppows(gmodulcp(concatsp(gmael(bnfz,8,5),gmael3(bnfz,8,4,2)),R),ell);
        !           858:   unitsellinv=grouppows(unitsell,-1);
        !           859:   unitsell=concatsp(unitsell,unitsellinv);
        !           860:   p1=unitsell;
        !           861:   for (i=2; i<=max((ell>>1),3); i++) p1=concatsp(p1,grouppows(unitsell,i));
        !           862:   unitsell=p1;
        !           863:   nmax=gnorml2(algtobasis(nfz,be));
        !           864:   lu=lg(unitsell)-1; ben=be;
        !           865:   do
        !           866:   {
        !           867:     be=ben;
        !           868:     for (i=1; i<=lu; i++)
        !           869:     {
        !           870:       p1=gmul(be,(GEN)unitsell[i]);
        !           871:       p2=gnorml2(algtobasis(nfz,p1));
        !           872:       if (gcmp(p2,nmax)<0)
        !           873:       {
        !           874:        nmax=p2; ben=p1;
        !           875:       }
        !           876:     }
        !           877:   }
        !           878:   while (!gegal(be,ben));
        !           879:   return be;
        !           880: }
        !           881:
        !           882: static GEN
        !           883: testx(GEN subgroup, GEN X, long prec)
        !           884: {
        !           885:   long i,v;
        !           886:   GEN be,polrelbe,p1;
        !           887:
        !           888: /* in alg. 5.3.18., C=nbcol */
        !           889:   X=gmul(unmodell,X);
        !           890:   if (gcmp0(X)) return gzero;
        !           891:   for (i=dv+1; i<=nbcol; i++) if (gcmp0((GEN)X[i])) return gzero;
        !           892:   for (i=1; i<=lSml2; i++)
        !           893:     if (gcmp0(gmul((GEN)vecMsup[i],X))) return gzero;
        !           894:   be=gun;
        !           895:   for (i=1; i<=nbcol; i++) be=gmul(be,gpui((GEN)vecw[i],lift((GEN)X[i]),0));
        !           896:   if (DEBUGLEVEL>=2) { fprintferr("reducing beta = "); outerr(be); }
        !           897:   be=reducebeta(be);
        !           898:   if (DEBUGLEVEL>=2) { fprintferr("beta reduced = "); outerr(be); }
        !           899:   polrelbe=computepolrelbeta(be);
        !           900:   v=varn(polrelbe);
        !           901:   p1=unifpol((GEN)bnf[7],polrelbe,0);
        !           902:   p1=denom(gtovec(p1));
        !           903:   polrelbe=gsubst(polrelbe,v,gdiv(polx[v],p1));
        !           904:   polrelbe=gmul(polrelbe,gpuigs(p1,degree(polrelbe)));
        !           905:   if (DEBUGLEVEL>=2) { fprintferr("polrelbe = "); outerr(polrelbe); }
        !           906:   p1=rnfconductor(bnf,polrelbe,prec);
        !           907:   if (!gegal((GEN)p1[1],module)) return gzero;
        !           908:   if (!gegal((GEN)p1[3],subgroup)) return gzero;
        !           909:   return polrelbe;
        !           910: }
        !           911:
        !           912: GEN
        !           913: rnfkummer(GEN bnr, GEN subgroup, long all, long prec)
        !           914: {
        !           915:   long i,j,av=avma,tetpil,llistpr,lfactell,e,vp,vd,ind;
        !           916:   GEN p1,p2,p3,p4,wk;
        !           917:   GEN rayclgp,bid,ideal;
        !           918:   GEN kk,clgp,fununits,torsunit,vecB,vecC,Tc,Tv,P;
        !           919:   GEN Q,Qt,idealz,gothf,factgothf,listpr,listex,factell,pp,p,pr,z,vecnul;
        !           920:   GEN M,al,K,Msup,X,finalresult,y,cycpro;
        !           921:
        !           922:   checkbnrgen(bnr);
        !           923:   wk=gmael4(bnr,1,8,4,1);
        !           924:   if (all) gell=stoi(all);
        !           925:   else
        !           926:   {
        !           927:     if (!gcmp0(subgroup)) gell=det(subgroup);
        !           928:     else gell=det(diagonal(gmael(bnr,5,2)));
        !           929:   }
        !           930:   if (gcmp1(gell)) { avma = av; return polx[varn(gmael3(bnr,1,7,1))]; }
        !           931:   if (!isprime(gell)) err(impl,"kummer for composite relative degree");
        !           932:   if (divise(wk,gell))
        !           933:     return gerepileupto(av,rnfkummersimple(bnr,subgroup,all,prec));
        !           934:   if (all && gcmp0(subgroup))
        !           935:     err(talker,"kummer when zeta not in K requires a specific subgroup");
        !           936:   bnf=(GEN)bnr[1];
        !           937:   nf=(GEN)bnf[7];
        !           938:   polnf=(GEN)nf[1]; vnf=varn(polnf); degK=degree(polnf);
        !           939:   if (!vnf) err(talker,"main variable in kummer must not be x");
        !           940:       /* step 7 */
        !           941:   p1=conductor(bnr,subgroup,2,prec);
        !           942:       /* fin step 7 */
        !           943:   bnr=(GEN)p1[2];
        !           944:   rayclgp=(GEN)bnr[5];
        !           945:   raycyc=(GEN)rayclgp[2]; lraycyc=lg(raycyc)-1;
        !           946:   bid=(GEN)bnr[2]; module=(GEN)bid[1];
        !           947:   ideal=(GEN)module[1];
        !           948:   if (gcmp0(subgroup)) subgroup = diagonal(raycyc);
        !           949:   ell=itos(gell);
        !           950:   unmodell=gmodulss(1,ell);
        !           951:       /* step 1 of alg 5.3.5. */
        !           952:   if (DEBUGLEVEL>=3) { fprintferr("Step 1\n"); flusherr(); }
        !           953:   phiell=cyclo(ell,vnf);
        !           954:   p1=(GEN)compositum2(polnf,phiell)[1];
        !           955:   kk=(GEN)p1[4];
        !           956:   R=(GEN)p1[1];
        !           957:   A1=(GEN)p1[2];
        !           958:   A2=(GEN)p1[3];
        !           959:   if (signe(leadingcoeff(R))<0)
        !           960:   {
        !           961:     R=gneg_i(R); A1=gmodulcp(lift(A1),R); A2=gmodulcp(lift(A2),R);
        !           962:   }
        !           963:       /* step 2 */
        !           964:   if (DEBUGLEVEL>=3) { fprintferr("Step 2\n"); flusherr(); }
        !           965:   degKz=degree(R);
        !           966:   m=degKz/degK;
        !           967:   d=(ell-1)/m;
        !           968:   g=lift(gpuigs(gener(gell),d));
        !           969:   if (gcmp1(lift(gpuigs(gmodulcp(g,stoi(ell*ell)),m)))) g=addsi(ell,g);
        !           970:       /* step reduction of R */
        !           971:   if (degKz<=20)
        !           972:   {
        !           973:     GEN A3,A3rev;
        !           974:
        !           975:     if (DEBUGLEVEL>=3) { fprintferr("Step reduction\n"); flusherr(); }
        !           976:     p1=polredabs2(R,prec);
        !           977:     if (DEBUGLEVEL>=3) { fprintferr("polredabs = ");outerr((GEN)p1[1]); }
        !           978:     R=(GEN)p1[1];
        !           979:     A3=(GEN)p1[2];
        !           980:     A1=gsubst(lift(A1),vnf,A3);
        !           981:     A2=gsubst(lift(A2),vnf,A3);
        !           982:     A3rev=polymodrecip(A3);
        !           983:     U=gadd(gpui(A2,g,0),gmul(kk,A1));
        !           984:     U=gsubst(lift(A3rev),vnf,U);
        !           985:   }
        !           986:   else U=gadd(gpui(A2,g,0),gmul(kk,A1));
        !           987:       /* step 3 */
        !           988:       /* on peut factoriser disc(R) avec th. 2.1.6. */
        !           989:   if (DEBUGLEVEL>=3) { fprintferr("Step 3\n"); flusherr(); }
        !           990:   bnfz=bnfinit0(R,0,NULL,prec); nfz=(GEN)bnfz[7];
        !           991:   clgp=gmael(bnfz,8,1);
        !           992:   cyc=(GEN)clgp[2]; gc=lg(cyc)-1; gencyc=(GEN)clgp[3];
        !           993:   for (j=1; j<=gc && divise((GEN)cyc[j],gell); j++);
        !           994:   rc=j-1;
        !           995:   vecalpha=cgetg(gc+1,t_VEC);
        !           996:   for (j=1; j<=gc; j++)
        !           997:   {
        !           998:     p1 = idealpow(nfz,(GEN)gencyc[j],(GEN)cyc[j]);
        !           999:     p1 = (GEN)isprincipalgenforce(bnfz, p1)[2];
        !          1000:     vecalpha[j] = (long)basistoalg(nfz, p1);
        !          1001:   }
        !          1002:   fununits = check_units(bnfz,"rnfkummer");
        !          1003:   torsunit = gmael3(bnfz,8,4,2);
        !          1004:   ru=(degKz>>1)-1;
        !          1005:   rv=rc+ru+1;
        !          1006:   vselmer=cgetg(rv+1,t_VEC);
        !          1007:   for (j=1; j<=rc; j++) vselmer[j]=vecalpha[j];
        !          1008:   for (   ; j<rv; j++) vselmer[j]=(long)gmodulcp((GEN)fununits[j-rc],R);
        !          1009:   vselmer[rv]=(long)gmodulcp((GEN)torsunit,R);
        !          1010:       /* computation of the uu(j) (see remark 5.2.15.) */
        !          1011:   uu=cgetg(gc+1,t_VEC);
        !          1012:   for (j=1; j<=rc; j++) uu[j]=zero;
        !          1013:   for (   ; j<=gc; j++) uu[j]=(long)lift(ginv(gmodulcp((GEN)cyc[j],gell)));
        !          1014:     /* step 4 */
        !          1015:   if (DEBUGLEVEL>=3) { fprintferr("Step 4\n"); flusherr(); }
        !          1016:   vecB=cgetg(rc+1,t_VEC);
        !          1017:   Tc=cgetg(rc+1,t_MAT);
        !          1018:   for (j=1; j<=rc; j++)
        !          1019:   {
        !          1020:     p1=isprincipalell(tauofideal((GEN)gencyc[j]));
        !          1021:     vecB[j]=p1[2];
        !          1022:     Tc[j]=p1[1];
        !          1023:   }
        !          1024:   p1=cgetg(m,t_VEC);
        !          1025:   p1[1]=(long)idmat(rc);
        !          1026:   for (j=2; j<=m-1; j++) p1[j]=lmul((GEN)p1[j-1],Tc);
        !          1027:   p2=cgetg(rc+1,t_VEC);
        !          1028:   for (j=1; j<=rc; j++) p2[j]=un;
        !          1029:   p3=vecB;
        !          1030:   for (j=1; j<=m-1; j++)
        !          1031:   {
        !          1032:     p3=gsubst(lift(p3),vnf,U);
        !          1033:     p4=groupproduct(grouppows(p3,(j*d)%ell),(GEN)p1[m-j]);
        !          1034:     for (i=1; i<=rc; i++) p2[i]=lmul((GEN)p2[i],(GEN)p4[i]);
        !          1035:   }
        !          1036:   vecC=p2;
        !          1037:       /* step 5 */
        !          1038:   if (DEBUGLEVEL>=3) { fprintferr("Step 5\n"); flusherr(); }
        !          1039:   Tv=cgetg(rv+1,t_MAT);
        !          1040:   for (j=1; j<=rv; j++)
        !          1041:   {
        !          1042:     Tv[j]=isvirtualunit(gsubst(lift((GEN)vselmer[j]),vnf,U))[1];
        !          1043:     if (DEBUGLEVEL>=3) { fprintferr("   %ld\n",j); flusherr(); }
        !          1044:   }
        !          1045:   P=lift(ker(gmul(unmodell,gsub(Tv,gmul(g,idmat(rv))))));
        !          1046:   dv=lg(P)-1;
        !          1047:   vecw=cgetg(dv+1,t_VEC);
        !          1048:   for (j=1; j<=dv; j++)
        !          1049:   {
        !          1050:     p1=gun;
        !          1051:     for (i=1; i<=rv; i++) p1=gmul(p1,gpui((GEN)vselmer[i],gcoeff(P,i,j),0));
        !          1052:     vecw[j]=(long)p1;
        !          1053:   }
        !          1054:       /* step 6 */
        !          1055:   if (DEBUGLEVEL>=3) { fprintferr("Step 6\n"); flusherr(); }
        !          1056:   Q=ker(gmul(unmodell,gsub(gtrans(Tc),gmul(g,idmat(rc)))));
        !          1057:   Qt=gtrans(Q);
        !          1058:   dc=lg(Q)-1;
        !          1059:       /* step 7 done above */
        !          1060:       /* step 8 */
        !          1061:   if (DEBUGLEVEL>=3) { fprintferr("Step 7 and 8\n"); flusherr(); }
        !          1062:   idealz=lifttokz(ideal);
        !          1063:   computematexpoteta1();
        !          1064:   if (!divise(idealnorm(nf,ideal),gell)) gothf=idealz;
        !          1065:   else
        !          1066:   {
        !          1067:     GEN subgroupz;
        !          1068:
        !          1069:     polrel=computepolrel();
        !          1070:     steinitzZk=steinitzaux(idmat(degKz));
        !          1071:     subgroupz=invimsubgroup(bnr,subgroup,idealz,prec);
        !          1072:     gothf=conductor(bnrz,subgroupz,0,prec);
        !          1073:   }
        !          1074:       /* step 9 */
        !          1075:   if (DEBUGLEVEL>=3) { fprintferr("Step 9\n"); flusherr(); }
        !          1076:   factgothf=idealfactor(nfz,gothf);
        !          1077:   listpr=(GEN)factgothf[1]; listex=(GEN)factgothf[2]; llistpr=lg(listpr)-1;
        !          1078:   factell=primedec(nfz,gell); lfactell=lg(factell)-1;
        !          1079:       /* step 10 and step 11 */
        !          1080:   if (DEBUGLEVEL>=3) { fprintferr("Step 10 and 11\n"); flusherr(); }
        !          1081:   Sm=cgetg(1,t_VEC);Sml1=cgetg(1,t_VEC);Sml2=cgetg(1,t_VEC);Sl=cgetg(1,t_VEC);
        !          1082:   ESml2=cgetg(1,t_VEC);
        !          1083:   for (i=1; i<=llistpr; i++)
        !          1084:   {
        !          1085:     pp=cgetg(2,t_VEC); pp[1]=listpr[i];
        !          1086:     p=gmael(pp,1,1); e=itos(gmael(pp,1,3));
        !          1087:     if (!gegal(p,gell))
        !          1088:     {
        !          1089:       if (!gcmp1((GEN)listex[i])) { avma=av; return gzero; }
        !          1090:       if (!isconjinprimelist(Sm,(GEN)pp[1])) Sm=concatsp(Sm,pp);
        !          1091:     }
        !          1092:     else
        !          1093:     {
        !          1094:       vp=itos((GEN)listex[i]); vd=(vp-1)*(ell-1)-ell*e;
        !          1095:       if (vd > 0) { avma=av; return gzero; }
        !          1096:       if (vd==0)
        !          1097:       {
        !          1098:        if (!isconjinprimelist(Sml1,(GEN)pp[1])) Sml1=concatsp(Sml1,pp);
        !          1099:       }
        !          1100:       else
        !          1101:       {
        !          1102:        if (vp==1) { avma=av; return gzero; }
        !          1103:         if (!isconjinprimelist(Sml2,(GEN)pp[1]))
        !          1104:         {
        !          1105:           Sml2=concatsp(Sml2,pp); ESml2=concatsp(ESml2,(GEN)listex[i]);
        !          1106:         }
        !          1107:       }
        !          1108:     }
        !          1109:   }
        !          1110:   for (i=1; i<=lfactell; i++)
        !          1111:     if (!idealval(nfz,gothf,(GEN)factell[i]))
        !          1112:     {
        !          1113:       pp=cgetg(2,t_VEC); pp[1]=factell[i];
        !          1114:       if (!isconjinprimelist(Sl,(GEN)pp[1])) Sl=concatsp(Sl,pp);
        !          1115:     }
        !          1116:   lSml2=lg(Sml2)-1; lSl=lg(Sl)-1; lSl2=lSml2+lSl;
        !          1117:   Sp=concatsp(Sm,Sml1); lSp=lg(Sp)-1;
        !          1118:       /* step 12 */
        !          1119:   if (DEBUGLEVEL>=3) { fprintferr("Step 12\n"); flusherr(); }
        !          1120:   vecbetap=cgetg(lSp+1,t_VEC);
        !          1121:   vecalphap=cgetg(lSp+1,t_VEC);
        !          1122:   matP=cgetg(lSp+1,t_MAT);
        !          1123:   gg=gmodulcp(g,gell);
        !          1124:   for (j=1; j<=lSp; j++)
        !          1125:   {
        !          1126:     p1=isprincipalell((GEN)Sp[j]);
        !          1127:     matP[j]=p1[1];
        !          1128:     p2=gun;
        !          1129:     for (i=1; i<=rc; i++)
        !          1130:       p2=gmul(p2,gpui((GEN)vecC[i],gmael(p1,1,i),0));
        !          1131:     p3=gdiv((GEN)p1[2],p2); vecbetap[j]=(long)p3;
        !          1132:     p2=gun;
        !          1133:     for (i=0; i<m; i++)
        !          1134:     {
        !          1135:       p2=gmul(p2,gpui(p3,lift(gpuigs(gg,m-1-i)),0));
        !          1136:       if (i<m-1) p3=gsubst(lift(p3),vnf,U);
        !          1137:     }
        !          1138:     vecalphap[j]=(long)p2;
        !          1139:   }
        !          1140:       /* step 13 */
        !          1141:   if (DEBUGLEVEL>=3) { fprintferr("Step 13\n"); flusherr(); }
        !          1142:   nbcol=lSp+dv;
        !          1143:   vecw=concatsp(vecw,vecbetap);
        !          1144:   listmod=cgetg(lSl2+1,t_VEC);
        !          1145:   listunif=cgetg(lSl2+1,t_VEC);
        !          1146:   listprSp=cgetg(lSl2+1,t_VEC);
        !          1147:   for (j=1; j<=lSml2; j++)
        !          1148:   {
        !          1149:     pr=(GEN)Sml2[j]; z=addsi(1,mulsi(ell,divis((GEN)pr[3],ell-1)));
        !          1150:     listmod[j]=(long)idealpow(nfz,pr,gsub(z,(GEN)ESml2[j]));
        !          1151:     listunif[j]=(long)basistoalg(nfz,gdiv((GEN)pr[5],(GEN)pr[1]));
        !          1152:     listprSp[j]=(long)pr;
        !          1153:   }
        !          1154:   for (   ; j<=lSl2; j++)
        !          1155:   {
        !          1156:     pr=(GEN)Sl[j-lSml2]; z=mulsi(ell,divis((GEN)pr[3],ell-1));
        !          1157:     listmod[j]=(long)idealpow(nfz,pr,z);
        !          1158:     listunif[j]=(long)basistoalg(nfz,gdiv((GEN)pr[5],(GEN)pr[1]));
        !          1159:     listprSp[j]=(long)pr;
        !          1160:   }
        !          1161:       /* step 14 and step 15 */
        !          1162:   if (DEBUGLEVEL>=3) { fprintferr("Step 14 and 15\n"); flusherr(); }
        !          1163:   listbid=cgetg(lSl2+1,t_VEC);
        !          1164:   listellrank=cgetg(lSl2+1,t_VEC);
        !          1165:   for (i=1; i<=lSl2; i++)
        !          1166:   {
        !          1167:     listbid[i]=(long)zidealstarinitgen(nfz,(GEN)listmod[i]);
        !          1168:     cycpro=gmael3(listbid,i,2,2);
        !          1169:     for (j=1; (j<lg(cycpro)) && (divise((GEN)cycpro[j],gell)); j++);
        !          1170:     listellrank[i]=lstoi(j-1);
        !          1171:   }
        !          1172:   mginv=lift(gmulsg(m,ginv(gg)));
        !          1173:   vecnul=cgetg(dc+1,t_COL); for (i=1; i<=dc; i++) vecnul[i]=zero;
        !          1174:   M=cgetg(nbcol+1,t_MAT);
        !          1175:   for (j=1; j<=dv; j++)
        !          1176:   {
        !          1177:     p1=cgetg(1,t_COL);
        !          1178:     al=(GEN)vecw[j];
        !          1179:     for (ind=1; ind<=lSl2; ind++) p1=concatsp(p1,ideallogaux(ind,al));
        !          1180:     p1=gmul(mginv,p1);
        !          1181:     M[j]=(long)concatsp(p1,vecnul);
        !          1182:   }
        !          1183:   for (   ; j<=nbcol; j++)
        !          1184:   {
        !          1185:     p1=cgetg(1,t_COL);
        !          1186:     al=(GEN)vecalphap[j-dv];
        !          1187:     for (ind=1; ind<=lSl2; ind++) p1=concatsp(p1,ideallogaux(ind,al));
        !          1188:     M[j]=(long)concatsp(p1,gmul(Qt,(GEN)matP[j-dv]));
        !          1189:   }
        !          1190:   M=gmul(unmodell,M);
        !          1191:       /* step 16 */
        !          1192:   if (DEBUGLEVEL>=3) { fprintferr("Step 16\n"); flusherr(); }
        !          1193:   K=ker(M);
        !          1194:   dK=lg(K)-1;
        !          1195:   if (!dK) { avma=av; return gzero; }
        !          1196:       /* step 17 */
        !          1197:   if (DEBUGLEVEL>=3) { fprintferr("Step 17\n"); flusherr(); }
        !          1198:   listmodsup=cgetg(lSml2+1,t_VEC);
        !          1199:   listbidsup=cgetg(lSml2+1,t_VEC);
        !          1200:   listellranksup=cgetg(lSml2+1,t_VEC);
        !          1201:   for (j=1; j<=lSml2; j++)
        !          1202:   {
        !          1203:     listmodsup[j]=(long)idealmul(nfz,(GEN)listprSp[j],(GEN)listmod[j]);
        !          1204:     listbidsup[j]=(long)zidealstarinitgen(nfz,(GEN)listmodsup[j]);
        !          1205:     cycpro=gmael3(listbidsup,j,2,2);
        !          1206:     for (i=1; (i<lg(cycpro)) && (divise((GEN)cycpro[i],gell)); i++);
        !          1207:     listellranksup[j]=lstoi(i-1);
        !          1208:   }
        !          1209:   vecMsup=cgetg(lSml2+1,t_VEC);
        !          1210:   for (i=1; i<=lSml2; i++)
        !          1211:   {
        !          1212:     Msup=cgetg(nbcol+1,t_MAT); vecMsup[i]=(long)Msup;
        !          1213:     for (j=1; j<=dv; j++) Msup[j]=lmul(mginv,ideallogauxsup(i,(GEN)vecw[j]));
        !          1214:     for (   ; j<=nbcol; j++) Msup[j]=(long)ideallogauxsup(i,(GEN)vecalphap[j-dv]);
        !          1215:   }
        !          1216:       /* step 18 */
        !          1217:   if (DEBUGLEVEL>=3) { fprintferr("Step 18\n"); flusherr(); }
        !          1218:   do
        !          1219:   {
        !          1220:     y=cgetg(dK,t_VEC);
        !          1221:     for (i=1; i<dK; i++) y[i]=zero;
        !          1222:        /* step 19 */
        !          1223:     LAB19:
        !          1224:     X=(GEN)K[dK];
        !          1225:     for (j=1; j<dK; j++) X=gadd(X,gmul((GEN)y[j],(GEN)K[j]));
        !          1226:     finalresult=testx(subgroup,X,prec);
        !          1227:     if (!gcmp0(finalresult))
        !          1228:     {
        !          1229:       tetpil=avma; return gerepile(av,tetpil,gcopy(finalresult));
        !          1230:     }
        !          1231:        /* step 20 */
        !          1232:     i=dK;
        !          1233:        /* step 21 */
        !          1234:     LAB21:
        !          1235:     i--;
        !          1236:        /* step 22 */
        !          1237:     if (i)
        !          1238:     {
        !          1239:       y[i]=(long)addsi(1,(GEN)y[i]);
        !          1240:       if (i<dK-1) y[i+1]=zero;
        !          1241:       if (cmpii((GEN)y[i],gell) >= 0) goto LAB21; else goto LAB19;
        !          1242:     }
        !          1243:     dK--;
        !          1244:   }
        !          1245:   while (dK);
        !          1246:   avma=av; return gzero;
        !          1247: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>