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

Annotation of OpenXM_contrib/pari-2.2/src/modules/kummer.c, Revision 1.1.1.1

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

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