[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

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>