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

Annotation of OpenXM_contrib/pari-2.2/src/basemath/base3.c, Revision 1.1

1.1     ! noro        1: /* $Id: base3.c,v 1.41 2001/10/01 12:11:29 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: /*                       BASIC NF OPERATIONS                       */
        !            19: /*                                                                 */
        !            20: /*******************************************************************/
        !            21: #include "pari.h"
        !            22: extern int ker_trivial_mod_p(GEN x, GEN p);
        !            23: extern long ideal_is_zk(GEN ideal,long N);
        !            24: extern GEN famat_ideallog(GEN nf, GEN g, GEN e, GEN bid);
        !            25: extern GEN famat_to_nf_modidele(GEN nf, GEN g, GEN e, GEN bid);
        !            26: extern GEN hnfall_i(GEN A, GEN *ptB, long remove);
        !            27: extern GEN vconcat(GEN A, GEN B);
        !            28:
        !            29: /*******************************************************************/
        !            30: /*                                                                 */
        !            31: /*                OPERATIONS OVER NUMBER FIELD ELEMENTS.           */
        !            32: /*     These are always represented as column vectors over the     */
        !            33: /*     integral basis nf[7]                                        */
        !            34: /*                                                                 */
        !            35: /*******************************************************************/
        !            36:
        !            37: int
        !            38: isnfscalar(GEN x)
        !            39: {
        !            40:   long lx=lg(x),i;
        !            41:
        !            42:   for (i=2; i<lx; i++)
        !            43:     if (!gcmp0((GEN)x[i])) return 0;
        !            44:   return 1;
        !            45: }
        !            46:
        !            47: static GEN
        !            48: scal_mul(GEN nf, GEN x, GEN y, long ty)
        !            49: {
        !            50:   long av=avma, tetpil;
        !            51:   GEN p1;
        !            52:
        !            53:   if (!is_extscalar_t(ty))
        !            54:   {
        !            55:     if (ty!=t_COL) err(typeer,"nfmul");
        !            56:     y = gmul((GEN)nf[7],y);
        !            57:   }
        !            58:   p1 = gmul(x,y); tetpil=avma;
        !            59:   return gerepile(av,tetpil,algtobasis(nf,p1));
        !            60: }
        !            61:
        !            62: static GEN
        !            63: get_tab(GEN nf, long *N)
        !            64: {
        !            65:   GEN tab = (typ(nf) == t_MAT)? nf: (GEN)nf[9];
        !            66:   *N = lg(tab[1])-1; return tab;
        !            67: }
        !            68:
        !            69: /* product of x and y in nf */
        !            70: GEN
        !            71: element_mul(GEN nf, GEN x, GEN y)
        !            72: {
        !            73:   long av,i,j,k,N,tx,ty;
        !            74:   GEN s,v,c,p1,tab;
        !            75:
        !            76:   if (x == y) return element_sqr(nf,x);
        !            77:
        !            78:   tx=typ(x); ty=typ(y);
        !            79:   nf=checknf(nf); tab = get_tab(nf, &N);
        !            80:   if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_mul");
        !            81:   if (ty==t_POLMOD) y=checknfelt_mod(nf,y,"element_mul");
        !            82:   if (is_extscalar_t(tx)) return scal_mul(nf,x,y,ty);
        !            83:   if (is_extscalar_t(ty)) return scal_mul(nf,y,x,tx);
        !            84:   if (isnfscalar(x)) return gmul((GEN)x[1],y);
        !            85:   if (isnfscalar(y)) return gmul((GEN)y[1],x);
        !            86:
        !            87:   v=cgetg(N+1,t_COL); av=avma;
        !            88:   for (k=1; k<=N; k++)
        !            89:   {
        !            90:     if (k == 1)
        !            91:       s = gmul((GEN)x[1],(GEN)y[1]);
        !            92:     else
        !            93:       s = gadd(gmul((GEN)x[1],(GEN)y[k]),
        !            94:                gmul((GEN)x[k],(GEN)y[1]));
        !            95:     for (i=2; i<=N; i++)
        !            96:     {
        !            97:       c=gcoeff(tab,k,(i-1)*N+i);
        !            98:       if (signe(c))
        !            99:       {
        !           100:         p1 = gmul((GEN)x[i],(GEN)y[i]);
        !           101:        if (!gcmp1(c)) p1 = gmul(p1,c);
        !           102:        s = gadd(s, p1);
        !           103:       }
        !           104:       for (j=i+1; j<=N; j++)
        !           105:       {
        !           106:        c=gcoeff(tab,k,(i-1)*N+j);
        !           107:        if (signe(c))
        !           108:        {
        !           109:           p1 = gadd(gmul((GEN)x[i],(GEN)y[j]),
        !           110:                     gmul((GEN)x[j],(GEN)y[i]));
        !           111:          if (!gcmp1(c)) p1 = gmul(p1,c);
        !           112:           s = gadd(s, p1);
        !           113:        }
        !           114:       }
        !           115:     }
        !           116:     v[k]=(long)gerepileupto(av,s); av=avma;
        !           117:   }
        !           118:   return v;
        !           119: }
        !           120:
        !           121: /* inverse of x in nf */
        !           122: GEN
        !           123: element_inv(GEN nf, GEN x)
        !           124: {
        !           125:   long av=avma,i,N,tx=typ(x);
        !           126:   GEN p1,p;
        !           127:
        !           128:   nf=checknf(nf); N=degpol(nf[1]);
        !           129:   if (is_extscalar_t(tx))
        !           130:   {
        !           131:     if (tx==t_POLMOD) checknfelt_mod(nf,x,"element_inv");
        !           132:     else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
        !           133:     return gerepileupto(av, algtobasis(nf, ginv(x)));
        !           134:   }
        !           135:   if (isnfscalar(x))
        !           136:   {
        !           137:     p1=cgetg(N+1,t_COL); p1[1]=linv((GEN)x[1]);
        !           138:     for (i=2; i<=N; i++) p1[i]=lcopy((GEN)x[i]);
        !           139:     return p1;
        !           140:   }
        !           141:   p = NULL;
        !           142:   for (i=1; i<=N; i++)
        !           143:     if (typ(x[i])==t_INTMOD)
        !           144:     {
        !           145:       p = gmael(x,i,1);
        !           146:       x = lift(x); break;
        !           147:     }
        !           148:   p1 = ginvmod(gmul((GEN)nf[7],x), (GEN)nf[1]);
        !           149:   p1 = algtobasis_intern(nf,p1);
        !           150:
        !           151:   if (p) p1 = FpV(p1, p);
        !           152:   return gerepileupto(av,p1);
        !           153: }
        !           154:
        !           155: /* quotient of x and y in nf */
        !           156: GEN
        !           157: element_div(GEN nf, GEN x, GEN y)
        !           158: {
        !           159:   long av=avma,i,N,tx=typ(x),ty=typ(y);
        !           160:   GEN p1,p;
        !           161:
        !           162:   nf=checknf(nf); N=degpol(nf[1]);
        !           163:   if (tx==t_POLMOD) checknfelt_mod(nf,x,"element_div");
        !           164:   else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
        !           165:
        !           166:   if (ty==t_POLMOD) checknfelt_mod(nf,y,"element_div");
        !           167:   else if (ty==t_POL) y=gmodulcp(y,(GEN)nf[1]);
        !           168:
        !           169:   if (is_extscalar_t(tx))
        !           170:   {
        !           171:     if (is_extscalar_t(ty)) p1=gdiv(x,y);
        !           172:     else
        !           173:     {
        !           174:       if (ty!=t_COL) err(typeer,"nfdiv");
        !           175:       p1=gdiv(x,gmodulcp(gmul((GEN)nf[7],y),(GEN)nf[1]));
        !           176:     }
        !           177:     return gerepileupto(av, algtobasis(nf,p1));
        !           178:   }
        !           179:   if (is_extscalar_t(ty))
        !           180:   {
        !           181:     if (tx!=t_COL) err(typeer,"nfdiv");
        !           182:     p1=gdiv(gmodulcp(gmul((GEN)nf[7],x),(GEN)nf[1]),y);
        !           183:     return gerepileupto(av, algtobasis(nf,p1));
        !           184:   }
        !           185:
        !           186:   if (isnfscalar(y)) return gdiv(x,(GEN)y[1]);
        !           187:   if (isnfscalar(x))
        !           188:   {
        !           189:     p1=element_inv(nf,y);
        !           190:     return gerepileupto(av, gmul((GEN)x[1],p1));
        !           191:   }
        !           192:
        !           193:   p = NULL;
        !           194:   for (i=1; i<=N; i++)
        !           195:     if (typ(x[i])==t_INTMOD)
        !           196:     {
        !           197:       p = gmael(x,i,1);
        !           198:       x = lift(x); break;
        !           199:     }
        !           200:   for (i=1; i<=N; i++)
        !           201:     if (typ(y[i])==t_INTMOD)
        !           202:     {
        !           203:       p1 = gmael(y,i,1);
        !           204:       if (p && !egalii(p,p1))
        !           205:         err(talker,"inconsistant prime moduli in element_inv");
        !           206:       y = lift(y); break;
        !           207:     }
        !           208:
        !           209:   p1 = gmul(gmul((GEN)nf[7],x), ginvmod(gmul((GEN)nf[7],y), (GEN)nf[1]));
        !           210:   p1 = algtobasis_intern(nf, gres(p1, (GEN)nf[1]));
        !           211:   if (p) p1 = FpV(p1,p);
        !           212:   return gerepileupto(av,p1);
        !           213: }
        !           214:
        !           215: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
        !           216: GEN
        !           217: element_muli(GEN nf, GEN x, GEN y)
        !           218: {
        !           219:   long av,i,j,k,N;
        !           220:   GEN p1,s,v,c,tab;
        !           221:
        !           222:   tab = get_tab(nf, &N);
        !           223:   v=cgetg(N+1,t_COL); av=avma;
        !           224:   for (k=1; k<=N; k++)
        !           225:   {
        !           226:     if (k == 1)
        !           227:       s = mulii((GEN)x[1],(GEN)y[1]);
        !           228:     else
        !           229:       s = addii(mulii((GEN)x[1],(GEN)y[k]),
        !           230:                 mulii((GEN)x[k],(GEN)y[1]));
        !           231:     for (i=2; i<=N; i++)
        !           232:     {
        !           233:       c=gcoeff(tab,k,(i-1)*N+i);
        !           234:       if (signe(c))
        !           235:       {
        !           236:         p1 = mulii((GEN)x[i],(GEN)y[i]);
        !           237:         if (!gcmp1(c)) p1 = mulii(p1,c);
        !           238:        s = addii(s,p1);
        !           239:       }
        !           240:       for (j=i+1; j<=N; j++)
        !           241:       {
        !           242:        c=gcoeff(tab,k,(i-1)*N+j);
        !           243:        if (signe(c))
        !           244:        {
        !           245:           p1 = addii(mulii((GEN)x[i],(GEN)y[j]),
        !           246:                      mulii((GEN)x[j],(GEN)y[i]));
        !           247:           if (!gcmp1(c)) p1 = mulii(p1,c);
        !           248:          s = addii(s,p1);
        !           249:        }
        !           250:       }
        !           251:     }
        !           252:     v[k]=(long) gerepileuptoint(av,s); av=avma;
        !           253:   }
        !           254:   return v;
        !           255: }
        !           256:
        !           257: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
        !           258: GEN
        !           259: element_sqri(GEN nf, GEN x)
        !           260: {
        !           261:   GEN p1,s,v,c,tab;
        !           262:   long av,i,j,k,N;
        !           263:
        !           264:   tab = get_tab(nf, &N);
        !           265:
        !           266:   v=cgetg(N+1,t_COL); av=avma;
        !           267:   for (k=1; k<=N; k++)
        !           268:   {
        !           269:     if (k == 1)
        !           270:       s = sqri((GEN)x[1]);
        !           271:     else
        !           272:       s = shifti(mulii((GEN)x[1],(GEN)x[k]), 1);
        !           273:     for (i=2; i<=N; i++)
        !           274:     {
        !           275:       c=gcoeff(tab,k,(i-1)*N+i);
        !           276:       if (signe(c))
        !           277:       {
        !           278:         p1 = sqri((GEN)x[i]);
        !           279:         if (!gcmp1(c)) p1 = mulii(p1,c);
        !           280:        s = addii(s,p1);
        !           281:       }
        !           282:       for (j=i+1; j<=N; j++)
        !           283:       {
        !           284:        c=gcoeff(tab,k,(i-1)*N+j);
        !           285:        if (signe(c))
        !           286:        {
        !           287:           p1 = shifti(mulii((GEN)x[i],(GEN)x[j]),1);
        !           288:           if (!gcmp1(c)) p1 = mulii(p1,c);
        !           289:          s = addii(s,p1);
        !           290:        }
        !           291:       }
        !           292:     }
        !           293:     v[k]=(long) gerepileuptoint(av,s); av=avma;
        !           294:   }
        !           295:   return v;
        !           296: }
        !           297:
        !           298: /* square of x in nf */
        !           299: GEN
        !           300: element_sqr(GEN nf, GEN x)
        !           301: {
        !           302:   long av = avma,i,j,k,N, tx = typ(x);
        !           303:   GEN p1,s,v,c,tab;
        !           304:
        !           305:   nf = checknf(nf);
        !           306:   tab = get_tab(nf, &N);
        !           307:   if (tx==t_POLMOD) x=checknfelt_mod(nf,x,"element_sqr");
        !           308:   if (is_extscalar_t(tx))
        !           309:     return gerepileupto(av, algtobasis(nf, gsqr(x)));
        !           310:   if (isnfscalar(x))
        !           311:   {
        !           312:     s=cgetg(N+1,t_COL); s[1]=lsqr((GEN)x[1]);
        !           313:     for (i=2; i<=N; i++) s[i]=lcopy((GEN)x[i]);
        !           314:     return s;
        !           315:   }
        !           316:   v=cgetg(N+1,t_COL); av = avma;
        !           317:   for (k=1; k<=N; k++)
        !           318:   {
        !           319:     if (k == 1)
        !           320:       s = gsqr((GEN)x[1]);
        !           321:     else
        !           322:       s = gmul2n(gmul((GEN)x[1],(GEN)x[k]), 1);
        !           323:     for (i=2; i<=N; i++)
        !           324:     {
        !           325:       c=gcoeff(tab,k,(i-1)*N+i);
        !           326:       if (signe(c))
        !           327:       {
        !           328:         p1 = gsqr((GEN)x[i]);
        !           329:        if (!gcmp1(c)) p1 = gmul(p1,c);
        !           330:         s = gadd(s,p1);
        !           331:       }
        !           332:       for (j=i+1; j<=N; j++)
        !           333:       {
        !           334:        c=gcoeff(tab,k,(i-1)*N+j);
        !           335:        if (signe(c))
        !           336:        {
        !           337:           p1 = gmul((GEN)x[i],(GEN)x[j]);
        !           338:          p1 = gcmp1(c)? gmul2n(p1,1): gmul(p1,shifti(c,1));
        !           339:          s = gadd(s,p1);
        !           340:        }
        !           341:       }
        !           342:     }
        !           343:     v[k]=(long)gerepileupto(av,s); av=avma;
        !           344:   }
        !           345:   return v;
        !           346: }
        !           347:
        !           348: /* Compute x^n in nf, left-shift binary powering */
        !           349: GEN
        !           350: element_pow(GEN nf, GEN x, GEN n)
        !           351: {
        !           352:   ulong av = avma;
        !           353:   long s,N,i,j,m;
        !           354:   GEN y,p1,cx;
        !           355:
        !           356:   if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
        !           357:   nf=checknf(nf); N=degpol(nf[1]);
        !           358:   s=signe(n); if (!s) return gscalcol_i(gun,N);
        !           359:   if (typ(x)!=t_COL) x=algtobasis(nf,x);
        !           360:
        !           361:   if (isnfscalar(x))
        !           362:   {
        !           363:     y = gscalcol_i(gun,N);
        !           364:     y[1] = (long)powgi((GEN)x[1],n); return y;
        !           365:   }
        !           366:   cx = content(x);
        !           367:   if (gcmp1(cx)) cx = NULL; else x = gdiv(x,cx);
        !           368:   p1 = n+2; m = *p1;
        !           369:   y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
        !           370:   for (i=lgefint(n)-2;;)
        !           371:   {
        !           372:     for (; j; m<<=1,j--)
        !           373:     {
        !           374:       y = element_sqr(nf, y);
        !           375:       if (m<0) y=element_mul(nf, y, x);
        !           376:     }
        !           377:     if (--i == 0) break;
        !           378:     m = *++p1, j = BITS_IN_LONG;
        !           379:   }
        !           380:   if (s<0) y = element_inv(nf, y);
        !           381:   if (cx) y = gmul(y, powgi(cx, n));
        !           382:   return av==avma? gcopy(y): gerepileupto(av,y);
        !           383: }
        !           384:
        !           385: /* x in Z^n, compute lift(x^n mod p) */
        !           386: GEN
        !           387: element_pow_mod_p(GEN nf, GEN x, GEN n, GEN p)
        !           388: {
        !           389:   ulong av = avma;
        !           390:   long s,N,i,j,m;
        !           391:   GEN y,p1;
        !           392:
        !           393:   if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
        !           394:   nf=checknf(nf); N=degpol(nf[1]);
        !           395:   s=signe(n); if (!s) return gscalcol_i(gun,N);
        !           396:   if (typ(x)!=t_COL) x=algtobasis(nf,x);
        !           397:
        !           398:   if (isnfscalar(x))
        !           399:   {
        !           400:     y = gscalcol_i(gun,N);
        !           401:     y[1] = (long)powmodulo((GEN)x[1],n,p); return y;
        !           402:   }
        !           403:   p1 = n+2; m = *p1;
        !           404:   y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
        !           405:   for (i=lgefint(n)-2;;)
        !           406:   {
        !           407:     for (; j; m<<=1,j--)
        !           408:     {
        !           409:       y = element_sqri(nf, y);
        !           410:       if (m<0) y=element_muli(nf, y, x);
        !           411:       y = FpV_red(y, p);
        !           412:     }
        !           413:     if (--i == 0) break;
        !           414:     m = *++p1, j = BITS_IN_LONG;
        !           415:   }
        !           416:   if (s<0)  y = FpV_red(element_inv(nf,y), p);
        !           417:   return av==avma? gcopy(y): gerepileupto(av,y);
        !           418: }
        !           419:
        !           420: /* x = I-th vector of the Z-basis of Z_K, in Z^n, compute lift(x^n mod p) */
        !           421: GEN
        !           422: element_powid_mod_p(GEN nf, long I, GEN n, GEN p)
        !           423: {
        !           424:   ulong av = avma;
        !           425:   long s,N,i,j,m;
        !           426:   GEN y,p1;
        !           427:
        !           428:   if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
        !           429:   nf=checknf(nf); N=degpol(nf[1]);
        !           430:   s=signe(n);
        !           431:   if (!s || I == 1) return gscalcol_i(gun,N);
        !           432:   p1 = n+2; m = *p1;
        !           433:   y = zerocol(N); y[I] = un;
        !           434:   j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
        !           435:   for (i=lgefint(n)-2;;)
        !           436:   {
        !           437:     for (; j; m<<=1,j--)
        !           438:     {
        !           439:       y = element_sqri(nf, y);
        !           440:       if (m<0) y=element_mulid(nf, y, I);
        !           441:       y = FpV_red(y, p);
        !           442:     }
        !           443:     if (--i == 0) break;
        !           444:     m = *++p1, j = BITS_IN_LONG;
        !           445:   }
        !           446:   if (s<0)  y = FpV_red(element_inv(nf,y), p);
        !           447:   return av==avma? gcopy(y): gerepileupto(av,y);
        !           448: }
        !           449:
        !           450: /* Outputs x.w_i, where w_i is the i-th elt of the integral basis */
        !           451: GEN
        !           452: element_mulid(GEN nf, GEN x, long i)
        !           453: {
        !           454:   long av,j,k,N;
        !           455:   GEN s,v,c,p1, tab;
        !           456:
        !           457:   if (i==1) return gcopy(x);
        !           458:   tab = get_tab(nf, &N);
        !           459:   if (typ(x) != t_COL || lg(x) != N+1) err(typeer,"element_mulid");
        !           460:   tab += (i-1)*N;
        !           461:   v=cgetg(N+1,t_COL); av=avma;
        !           462:   for (k=1; k<=N; k++)
        !           463:   {
        !           464:     s = gzero;
        !           465:     for (j=1; j<=N; j++)
        !           466:       if (signe(c = gcoeff(tab,k,j)) && !gcmp0(p1 = (GEN)x[j]))
        !           467:       {
        !           468:         if (!gcmp1(c)) p1 = gmul(p1,c);
        !           469:        s = gadd(s,p1);
        !           470:       }
        !           471:     v[k]=lpileupto(av,s); av=avma;
        !           472:   }
        !           473:   return v;
        !           474: }
        !           475:
        !           476: /* valuation of integer x, with resp. to prime ideal P above p.
        !           477:  * p.P^(-1) = b Z_K, v >= val_p(norm(x)), and N = deg(nf)
        !           478:  * [b may be given as the 'multiplication by b' matrix]
        !           479:  */
        !           480: long
        !           481: int_elt_val(GEN nf, GEN x, GEN p, GEN b, GEN *newx, long v)
        !           482: {
        !           483:   long i,k,w, N = degpol(nf[1]);
        !           484:   GEN r,a,y,mul;
        !           485:
        !           486:   if (typ(b) == t_MAT) mul = b;
        !           487:   else
        !           488:   {
        !           489:     mul = cgetg(N+1,t_MAT);
        !           490:     for (i=1; i<=N; i++) mul[i] = (long)element_mulid(nf,b,i);
        !           491:   }
        !           492:   y = cgetg(N+1, t_COL); /* will hold the new x */
        !           493:   x = dummycopy(x);
        !           494:   for(w=0; w<=v; w++)
        !           495:   {
        !           496:     for (i=1; i<=N; i++)
        !           497:     { /* compute (x.b)_i */
        !           498:       a = mulii((GEN)x[1], gcoeff(mul,i,1));
        !           499:       for (k=2; k<=N; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
        !           500:       /* is it divisible by p ? */
        !           501:       y[i] = ldvmdii(a,p,&r);
        !           502:       if (signe(r)) goto END;
        !           503:     }
        !           504:     r=x; x=y; y=r;
        !           505:   }
        !           506: END:
        !           507:   if (newx) *newx = x;
        !           508:   return w;
        !           509: }
        !           510:
        !           511: long
        !           512: element_val(GEN nf, GEN x, GEN vp)
        !           513: {
        !           514:   long av = avma,N,w,vcx,e;
        !           515:   GEN cx,p;
        !           516:
        !           517:   if (gcmp0(x)) return VERYBIGINT;
        !           518:   nf=checknf(nf); N=degpol(nf[1]);
        !           519:   checkprimeid(vp); p=(GEN)vp[1]; e=itos((GEN)vp[3]);
        !           520:   switch(typ(x))
        !           521:   {
        !           522:     case t_INT: case t_FRAC: case t_FRACN:
        !           523:       return ggval(x,p)*e;
        !           524:     case t_POLMOD: x = (GEN)x[2]; /* fall through */
        !           525:     case t_POL:
        !           526:       x = algtobasis_intern(nf,x); break;
        !           527:     case t_COL:
        !           528:       if (lg(x)==N+1) break;
        !           529:     default: err(typeer,"element_val");
        !           530:   }
        !           531:   if (isnfscalar(x)) return ggval((GEN)x[1],p)*e;
        !           532:
        !           533:   cx = content(x);
        !           534:   if (gcmp1(cx)) vcx=0; else { x = gdiv(x,cx); vcx = ggval(cx,p); }
        !           535:   w = int_elt_val(nf,x,p,(GEN)vp[5],NULL,VERYBIGINT);
        !           536:   avma=av; return w + vcx*e;
        !           537: }
        !           538:
        !           539: /* d = a multiple of norm(x), x integral */
        !           540: long
        !           541: element_val2(GEN nf, GEN x, GEN d, GEN vp)
        !           542: {
        !           543:   GEN p = (GEN)vp[1];
        !           544:   long av, v = ggval(d,p);
        !           545:
        !           546:   if (!v) return 0;
        !           547:   av=avma;
        !           548:   v = int_elt_val(nf,x,p,(GEN)vp[5],NULL,v);
        !           549:   avma=av; return v;
        !           550: }
        !           551:
        !           552: /* polegal without comparing variables */
        !           553: long
        !           554: polegal_spec(GEN x, GEN y)
        !           555: {
        !           556:   long i = lgef(x);
        !           557:
        !           558:   if (i != lgef(y)) return 0;
        !           559:   for (i--; i > 1; i--)
        !           560:     if (!gegal((GEN)x[i],(GEN)y[i])) return 0;
        !           561:   return 1;
        !           562: }
        !           563:
        !           564: GEN
        !           565: basistoalg(GEN nf, GEN x)
        !           566: {
        !           567:   long tx=typ(x),lx=lg(x),i;
        !           568:   GEN z;
        !           569:
        !           570:   nf=checknf(nf);
        !           571:   switch(tx)
        !           572:   {
        !           573:     case t_COL:
        !           574:       for (i=1; i<lx; i++)
        !           575:       {
        !           576:         long t = typ(x[i]);
        !           577:        if (is_matvec_t(t)) break;
        !           578:       }
        !           579:       if (i==lx)
        !           580:       {
        !           581:         z = cgetg(3,t_POLMOD);
        !           582:         z[1] = lcopy((GEN)nf[1]);
        !           583:        z[2] = lmul((GEN)nf[7],x); return z;
        !           584:       }
        !           585:       /* fall through */
        !           586:
        !           587:     case t_VEC: case t_MAT: z=cgetg(lx,tx);
        !           588:       for (i=1; i<lx; i++) z[i]=(long)basistoalg(nf,(GEN)x[i]);
        !           589:       return z;
        !           590:
        !           591:     case t_POLMOD:
        !           592:       if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
        !           593:        err(talker,"not the same number field in basistoalg");
        !           594:       return gcopy(x);
        !           595:     default: z=cgetg(3,t_POLMOD);
        !           596:       z[1]=lcopy((GEN)nf[1]);
        !           597:       z[2]=lmul(x,polun[varn(nf[1])]); return z;
        !           598:   }
        !           599: }
        !           600:
        !           601: /* return the (N-dimensional) vector of coeffs of p */
        !           602: GEN
        !           603: pol_to_vec(GEN x, long N)
        !           604: {
        !           605:   long i, l;
        !           606:   GEN z = cgetg(N+1,t_COL);
        !           607:   if (typ(x) != t_POL)
        !           608:   {
        !           609:     z[1] = (long)x;
        !           610:     for (i=2; i<=N; i++) z[i]=zero;
        !           611:     return z;
        !           612:   }
        !           613:   l = lgef(x)-1; x++;
        !           614:   for (i=1; i<l ; i++) z[i]=x[i];
        !           615:   for (   ; i<=N; i++) z[i]=zero;
        !           616:   return z;
        !           617: }
        !           618:
        !           619: /* gmul(A, pol_to_vec(x)), A matrix of compatible dimensions */
        !           620: GEN
        !           621: mulmat_pol(GEN A, GEN x)
        !           622: {
        !           623:   long i,l;
        !           624:   GEN z;
        !           625:   if (typ(x) != t_POL) return gscalcol(x, lg(A[1])-1);
        !           626:   l=lgef(x)-1; if (l == 1) return zerocol(lg(A[1])-1);
        !           627:   x++; z = gmul((GEN)x[1], (GEN)A[1]);
        !           628:   for (i=2; i<l ; i++)
        !           629:     if (!gcmp0((GEN)x[i])) z = gadd(z, gmul((GEN)x[i], (GEN)A[i]));
        !           630:   return z;
        !           631: }
        !           632:
        !           633: /* valid for scalars and polynomial, degree less than N.
        !           634:  * No garbage collecting. No check (SEGV for vectors).
        !           635:  */
        !           636: GEN
        !           637: algtobasis_intern(GEN nf, GEN x)
        !           638: {
        !           639:   GEN P = (GEN)nf[1];
        !           640:   long tx = typ(x), N = degpol(P);
        !           641:
        !           642:   if (tx==t_POLMOD) { x=(GEN)x[2]; tx=typ(x); }
        !           643:   if (tx==t_POL)
        !           644:   {
        !           645:     if (varn(x) != varn(P))
        !           646:       err(talker,"incompatible variables in algtobasis");
        !           647:     if (degpol(x) >= N) x = gres(x,P);
        !           648:     return mulmat_pol((GEN)nf[8], x);
        !           649:   }
        !           650:   return gscalcol(x,N);
        !           651: }
        !           652:
        !           653: GEN
        !           654: algtobasis(GEN nf, GEN x)
        !           655: {
        !           656:   long tx=typ(x),lx=lg(x),av=avma,i,N;
        !           657:   GEN z;
        !           658:
        !           659:   nf=checknf(nf);
        !           660:   switch(tx)
        !           661:   {
        !           662:     case t_VEC: case t_COL: case t_MAT:
        !           663:       z=cgetg(lx,tx);
        !           664:       for (i=1; i<lx; i++) z[i]=(long)algtobasis(nf,(GEN)x[i]);
        !           665:       return z;
        !           666:     case t_POLMOD:
        !           667:       if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
        !           668:        err(talker,"not the same number field in algtobasis");
        !           669:       x = (GEN)x[2]; /* fall through */
        !           670:     case t_POL:
        !           671:       return gerepileupto(av,algtobasis_intern(nf,x));
        !           672:
        !           673:     default: N=degpol(nf[1]); return gscalcol(x,N);
        !           674:   }
        !           675: }
        !           676:
        !           677: /* Given a and b in nf, gives an algebraic integer y in nf such that a-b.y
        !           678:  * is "small"
        !           679:  */
        !           680: GEN
        !           681: nfdiveuc(GEN nf, GEN a, GEN b)
        !           682: {
        !           683:   long av=avma, tetpil;
        !           684:   a = element_div(nf,a,b); tetpil=avma;
        !           685:   return gerepile(av,tetpil,ground(a));
        !           686: }
        !           687:
        !           688: /* Given a and b in nf, gives a "small" algebraic integer r in nf
        !           689:  * of the form a-b.y
        !           690:  */
        !           691: GEN
        !           692: nfmod(GEN nf, GEN a, GEN b)
        !           693: {
        !           694:   long av=avma,tetpil;
        !           695:   GEN p1=gneg_i(element_mul(nf,b,ground(element_div(nf,a,b))));
        !           696:   tetpil=avma; return gerepile(av,tetpil,gadd(a,p1));
        !           697: }
        !           698:
        !           699: /* Given a and b in nf, gives a two-component vector [y,r] in nf such
        !           700:  * that r=a-b.y is "small".
        !           701:  */
        !           702: GEN
        !           703: nfdivres(GEN nf, GEN a, GEN b)
        !           704: {
        !           705:   long av=avma,tetpil;
        !           706:   GEN p1,z, y = ground(element_div(nf,a,b));
        !           707:
        !           708:   p1=gneg_i(element_mul(nf,b,y)); tetpil=avma;
        !           709:   z=cgetg(3,t_VEC); z[1]=lcopy(y); z[2]=ladd(a,p1);
        !           710:   return gerepile(av,tetpil,z);
        !           711: }
        !           712:
        !           713: /*************************************************************************/
        !           714: /**                                                                    **/
        !           715: /**                          (Z_K/I)^*                                 **/
        !           716: /**                                                                    **/
        !           717: /*************************************************************************/
        !           718:
        !           719: /* return (column) vector of R1 signatures of x (coeff modulo 2)
        !           720:  * if arch = NULL, assume arch = [0,..0]
        !           721:  */
        !           722: GEN
        !           723: zsigne(GEN nf,GEN x,GEN arch)
        !           724: {
        !           725:   GEN _0,_1, vecsign, rac = (GEN)nf[6];
        !           726:   long i,j,l,e,B;
        !           727:   ulong av0, av;
        !           728:
        !           729:   if (!arch) return cgetg(1,t_COL);
        !           730:   av0 = avma;
        !           731:   l = lg(arch); vecsign = cgetg(l,t_COL);
        !           732:   _0 = gmodulss(0,2);
        !           733:   _1 = gmodulss(1,2); av = avma;
        !           734:   switch(typ(x))
        !           735:   {
        !           736:     case t_MAT: /* factorisation */
        !           737:     {
        !           738:       GEN t = (GEN)x[1], e = (GEN)x[2];
        !           739:       for (j=i=1; i<l; i++)
        !           740:         if (signe(arch[i])) vecsign[j++] = (long)_0;
        !           741:       setlg(vecsign,j);
        !           742:       if (j==1) { avma = av0; return cgetg(1, t_COL); }
        !           743:       for (i=1; i<lg(t); i++)
        !           744:       {
        !           745:         GEN p1 = (GEN)e[i];
        !           746:         if (mpodd(p1))
        !           747:           vecsign = gadd(vecsign, zsigne(nf,(GEN)t[i],arch));
        !           748:       }
        !           749:       return gerepileupto(av0, vecsign);
        !           750:     }
        !           751:     case t_COL: x = gmul((GEN)nf[7],x); break;
        !           752:     case t_POLMOD: x = (GEN)x[2];
        !           753:   }
        !           754:   if (gcmp0(x)) err(talker,"zero element in zsigne");
        !           755:
        !           756:   B = bit_accuracy(precision((GEN)rac[1]));
        !           757:   e = gexpo(x);
        !           758:   for (j=i=1; i<l; i++)
        !           759:     if (signe(arch[i]))
        !           760:     {
        !           761:       GEN y = poleval(x,(GEN)rac[i]);
        !           762:       if (e + gexpo((GEN)rac[i]) - gexpo(y) > B)
        !           763:         err(precer, "zsigne");
        !           764:       vecsign[j++] = (gsigne(y) > 0)? (long)_0: (long)_1;
        !           765:     }
        !           766:   avma = av; setlg(vecsign,j); return vecsign;
        !           767: }
        !           768:
        !           769: GEN
        !           770: lllreducemodmatrix(GEN x,GEN y)
        !           771: {
        !           772:   ulong av = avma;
        !           773:   GEN z = gmul(y,lllint(y));
        !           774:   return gerepileupto(av, reducemodinvertible(x, z));
        !           775: }
        !           776:
        !           777: /* for internal use...reduce x modulo (invertible) y */
        !           778: GEN
        !           779: reducemodinvertible(GEN x, GEN y)
        !           780: {
        !           781:   return gadd(x, gneg(gmul(y, ground(gauss(y,x)))));
        !           782: }
        !           783:
        !           784: /* Reduce column x modulo y in HNF */
        !           785: GEN
        !           786: colreducemodmat(GEN x, GEN y, GEN *Q)
        !           787: {
        !           788:   long i, l = lg(x);
        !           789:   GEN q;
        !           790:
        !           791:   if (Q) *Q = cgetg(l,t_COL);
        !           792:   if (l == 1) return cgetg(1,t_COL);
        !           793:   for (i = l-1; i>0; i--)
        !           794:   {
        !           795:     q = negi(gdivround((GEN)x[i], gcoeff(y,i,i)));
        !           796:     if (Q) (*Q)[i] = (long)q;
        !           797:     if (signe(q)) x = gadd(x, gmul(q, (GEN)y[i]));
        !           798:   }
        !           799:   return x;
        !           800: }
        !           801:
        !           802: /* for internal use...Reduce matrix x modulo matrix y */
        !           803: GEN
        !           804: reducemodmatrix(GEN x, GEN y)
        !           805: {
        !           806:   return reducemodHNF(x, hnfmod(y,detint(y)), NULL);
        !           807: }
        !           808:
        !           809: /* x = y Q + R */
        !           810: GEN
        !           811: reducemodHNF(GEN x, GEN y, GEN *Q)
        !           812: {
        !           813:   long lx = lg(x), i;
        !           814:   GEN R = cgetg(lx, t_MAT);
        !           815:   if (Q)
        !           816:   {
        !           817:     GEN q = cgetg(lx, t_MAT); *Q = q;
        !           818:     for (i=1; i<lx; i++) R[i] = (long)colreducemodmat((GEN)x[i],y,(GEN*)(q+i));
        !           819:   }
        !           820:   else
        !           821:     for (i=1; i<lx; i++)
        !           822:     {
        !           823:       long av = avma;
        !           824:       R[i] = lpileupto(av, colreducemodmat((GEN)x[i],y,NULL));
        !           825:     }
        !           826:   return R;
        !           827: }
        !           828:
        !           829: /* For internal use. Reduce x modulo ideal. We want a non-zero result */
        !           830: GEN
        !           831: nfreducemodideal(GEN nf,GEN x0,GEN ideal)
        !           832: {
        !           833:   GEN y = idealhermite(nf,ideal);
        !           834:   GEN x = colreducemodmat(x0, y, NULL);
        !           835:   if (gcmp0(x)) return (GEN)y[1];
        !           836:   return x == x0? gcopy(x) : x;
        !           837: }
        !           838:
        !           839: /* multiply y by t = 1 mod^* f such that sign(x) = sign(y) at arch = idele[2].
        !           840:  * If x == NULL, make y >> 0 at sarch */
        !           841: GEN
        !           842: set_sign_mod_idele(GEN nf, GEN x, GEN y, GEN idele, GEN sarch)
        !           843: {
        !           844:   GEN s,arch,gen;
        !           845:   long nba,i;
        !           846:   if (!sarch) return y;
        !           847:   gen = (GEN)sarch[2]; nba = lg(gen);
        !           848:   if (nba == 1) return y;
        !           849:
        !           850:   arch = (GEN)idele[2];
        !           851:   s = zsigne(nf,y,arch);
        !           852:   if (x) s = gadd(s, zsigne(nf,x,arch));
        !           853:   s = lift_intern(gmul((GEN)sarch[3],s));
        !           854:   for (i=1; i<nba; i++)
        !           855:     if (signe(s[i])) y = element_mul(nf,y,(GEN)gen[i]);
        !           856:   return y;
        !           857: }
        !           858:
        !           859: /* compute elt = x mod idele, with sign(elt) = sign(x) at arch */
        !           860: GEN
        !           861: nfreducemodidele(GEN nf,GEN x,GEN idele,GEN sarch)
        !           862: {
        !           863:   GEN y;
        !           864:
        !           865:   if (gcmp0(x)) return gcopy(x);
        !           866:   if (!sarch || typ(idele)!=t_VEC || lg(idele)!=3)
        !           867:     return nfreducemodideal(nf,x,idele);
        !           868:
        !           869:   y = nfreducemodideal(nf,x,(GEN)idele[1]);
        !           870:   y = set_sign_mod_idele(nf, x, y, idele, sarch);
        !           871:   return (gexpo(y) > gexpo(x))? x: y;
        !           872: }
        !           873:
        !           874: GEN
        !           875: element_powmodideal(GEN nf,GEN x,GEN k,GEN ideal)
        !           876: {
        !           877:   GEN y = gscalcol_i(gun,degpol(nf[1]));
        !           878:   for(;;)
        !           879:   {
        !           880:     if (mpodd(k)) y=element_mulmodideal(nf,x,y,ideal);
        !           881:     k=shifti(k,-1); if (!signe(k)) return y;
        !           882:     x = element_sqrmodideal(nf,x,ideal);
        !           883:   }
        !           884: }
        !           885:
        !           886: GEN
        !           887: element_powmodidele(GEN nf,GEN x,GEN k,GEN idele,GEN sarch)
        !           888: {
        !           889:   GEN y = element_powmodideal(nf,x,k,idele);
        !           890:   if (mpodd(k))
        !           891:   { if (!gcmp1(k)) y = set_sign_mod_idele(nf, x, y, idele, sarch); }
        !           892:   else
        !           893:   { if (!gcmp0(k)) y = set_sign_mod_idele(nf, NULL, y, idele, sarch); }
        !           894:   return y;
        !           895: }
        !           896:
        !           897: /* H relation matrix among row of generators g.  Let URV = D its SNF, newU R
        !           898:  * newV = newD its clean SNF (no 1 in Dnew). Return the diagonal of newD,
        !           899:  * newU and newUi such that  1/U = (newUi, ?).
        !           900:  * Rationale: let (G,0) = g Ui be the new generators then
        !           901:  * 0 = G U R --> G D = 0,  g = G newU  and  G = g newUi */
        !           902: GEN
        !           903: smithrel(GEN H, GEN *newU, GEN *newUi)
        !           904: {
        !           905:   GEN z,U,Ui,D,p1;
        !           906:   long l,c;
        !           907:
        !           908:   z = smith2(H);
        !           909:   U = (GEN)z[1];
        !           910:   D = (GEN)z[3]; l = lg(D);
        !           911:   for (c=1; c<l; c++)
        !           912:   {
        !           913:     p1 = gcoeff(D,c,c);
        !           914:     if (is_pm1(p1)) break;
        !           915:   }
        !           916:   if (newU) *newU = rowextract_i(U, 1, c-1);
        !           917:   if (newUi) {
        !           918:     Ui = ginv(U); setlg(Ui, c);
        !           919:     *newUi = reducemodHNF(Ui, H, NULL);
        !           920:   }
        !           921:   setlg(D, c); return mattodiagonal_i(D);
        !           922: }
        !           923:
        !           924: /* given 2 integral ideals x, y in HNF s.t x | y | x^2, compute the quotient
        !           925:    (1+x)/(1+y) in the form [[cyc],[gen],ux^-1]. */
        !           926: static GEN
        !           927: zidealij(GEN x, GEN y)
        !           928: {
        !           929:   GEN U,R,G,p1,cyc,z,xi;
        !           930:   long j,N;
        !           931:
        !           932:   xi = ginv(x); R = gmul(xi, y); /* relations between the 1 + x_i (HNF) */
        !           933:   cyc = smithrel(R, &U, &G);
        !           934:   N = lg(cyc); G = gmul(x,G); settyp(G, t_VEC); /* new generators */
        !           935:   for (j=1; j<N; j++)
        !           936:   {
        !           937:     p1 = (GEN)G[j];
        !           938:     p1[1] = laddsi(1, (GEN)p1[1]); /* 1 + g_j */
        !           939:   }
        !           940:   z = cgetg(4,t_VEC);
        !           941:   z[1] = (long)cyc;
        !           942:   z[2] = (long)G;
        !           943:   z[3] = lmul(U,xi); return z;
        !           944: }
        !           945:
        !           946: static GEN
        !           947: get_multab(GEN nf, GEN x)
        !           948: {
        !           949:   long lx = lg(x), i;
        !           950:   GEN multab = cgetg(lx, t_MAT);
        !           951:   for (i=1; i<lx; i++)
        !           952:     multab[i]=(long)element_mulid(nf,x,i);
        !           953:   return multab;
        !           954: }
        !           955:
        !           956: /* return mat * y mod prh */
        !           957: static GEN
        !           958: mul_matvec_mod_pr(GEN mat, GEN y, GEN prh)
        !           959: {
        !           960:   long av,i,j, lx = lg(mat);
        !           961:   GEN v, res = cgetg(lx,t_COL), p = gcoeff(prh,1,1);
        !           962:
        !           963:   av = avma; (void)new_chunk((lx-1) * lgefint(p));
        !           964:   v = zerocol(lx-1);
        !           965:   for (i=lx-1; i; i--)
        !           966:   {
        !           967:     GEN p1 = (GEN)v[i], t = (GEN)prh[i];
        !           968:     for (j=1; j<lx; j++)
        !           969:       p1 = addii(p1, mulii(gcoeff(mat,i,j), (GEN)y[j]));
        !           970:     p1 = modii(p1, p);
        !           971:     if (p1 != gzero && is_pm1(t[i]))
        !           972:     {
        !           973:       for (j=1; j<i; j++)
        !           974:         v[j] = lsubii((GEN)v[j], mulii(p1, (GEN)t[j]));
        !           975:       p1 = gzero;
        !           976:     }
        !           977:     if (p1 == gzero) /* intended */
        !           978:       res[i] = zero;
        !           979:     else
        !           980:       av = res[i] = (long)icopy_av(p1, (GEN)av);
        !           981:   }
        !           982:   avma = av; return res;
        !           983: }
        !           984:
        !           985: /* smallest integer n such that g0^n=x modulo p prime. Assume g0 has order q
        !           986:  * (p-1 if q = NULL) */
        !           987: GEN
        !           988: Fp_shanks(GEN x,GEN g0,GEN p, GEN q)
        !           989: {
        !           990:   long av=avma,av1,lim,lbaby,i,k,c;
        !           991:   GEN p1,smalltable,giant,perm,v,g0inv;
        !           992:
        !           993:   x = modii(x,p);
        !           994:   if (is_pm1(x) || egalii(p,gdeux)) { avma = av; return gzero; }
        !           995:   p1 = addsi(-1, p); if (!q) q = p1;
        !           996:   if (egalii(p1,x)) { avma = av; return shifti(q,-1); }
        !           997:   p1 = racine(q);
        !           998:   if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in Fp_shanks");
        !           999:   lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
        !          1000:   g0inv = mpinvmod(g0, p); p1 = x;
        !          1001:
        !          1002:   c = 3 * lgefint(p);
        !          1003:   for (i=1;;i++)
        !          1004:   {
        !          1005:     av1 = avma;
        !          1006:     if (is_pm1(p1)) { avma=av; return stoi(i-1); }
        !          1007:     smalltable[i]=(long)p1; if (i==lbaby) break;
        !          1008:     new_chunk(c); p1 = mulii(p1,g0inv);
        !          1009:     avma = av1; p1 = resii(p1, p);
        !          1010:   }
        !          1011:   giant = resii(mulii(x, mpinvmod(p1,p)), p);
        !          1012:   p1=cgetg(lbaby+1,t_VEC);
        !          1013:   perm = gen_sort(smalltable, cmp_IND | cmp_C, cmpii);
        !          1014:   for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
        !          1015:   smalltable=p1; p1=giant;
        !          1016:
        !          1017:   av1 = avma; lim=stack_lim(av1,2);
        !          1018:   for (k=1;;k++)
        !          1019:   {
        !          1020:     i=tablesearch(smalltable,p1,cmpii);
        !          1021:     if (i)
        !          1022:     {
        !          1023:       v=addis(mulss(lbaby-1,k),perm[i]);
        !          1024:       return gerepileuptoint(av,addsi(-1,v));
        !          1025:     }
        !          1026:     p1 = resii(mulii(p1,giant), p);
        !          1027:
        !          1028:     if (low_stack(lim, stack_lim(av1,2)))
        !          1029:     {
        !          1030:       if(DEBUGMEM>1) err(warnmem,"Fp_shanks, k = %ld", k);
        !          1031:       p1 = gerepileuptoint(av1, p1);
        !          1032:     }
        !          1033:   }
        !          1034: }
        !          1035:
        !          1036: /* Pohlig-Hellman */
        !          1037: GEN
        !          1038: Fp_PHlog(GEN a, GEN g, GEN p, GEN ord)
        !          1039: {
        !          1040:   ulong av = avma;
        !          1041:   GEN v,t0,a0,b,q,g_q,n_q,ginv0,qj,ginv;
        !          1042:   GEN fa, ex;
        !          1043:   long e,i,j,l;
        !          1044:
        !          1045:   if (!ord) ord = subis(p,1);
        !          1046:   if (typ(ord) == t_MAT)
        !          1047:   {
        !          1048:     fa = ord;
        !          1049:     ord= factorback(fa,NULL);
        !          1050:   }
        !          1051:   else
        !          1052:     fa = decomp(ord);
        !          1053:   if (typ(g) == t_INTMOD) g = lift_intern(g);
        !          1054:   ex = (GEN)fa[2];
        !          1055:   fa = (GEN)fa[1];
        !          1056:   l = lg(fa);
        !          1057:   ginv = mpinvmod(g,p);
        !          1058:   v = cgetg(l, t_VEC);
        !          1059:   for (i=1; i<l; i++)
        !          1060:   {
        !          1061:     q = (GEN)fa[i];
        !          1062:     e = itos((GEN)ex[i]);
        !          1063:     if (DEBUGLEVEL>5)
        !          1064:       fprintferr("Pohlig-Hellman: DL mod %Z^%ld\n",q,e);
        !          1065:     qj = new_chunk(e+1); qj[0] = un;
        !          1066:     for (j=1; j<=e; j++) qj[j] = lmulii((GEN)qj[j-1], q);
        !          1067:     t0 = divii(ord, (GEN)qj[e]);
        !          1068:     a0 = powmodulo(a, t0, p);
        !          1069:     ginv0 = powmodulo(ginv, t0, p); /* order q^e */
        !          1070:     g_q = powmodulo(g, divii(ord,q), p); /* order q */
        !          1071:     n_q = gzero;
        !          1072:     for (j=0; j<e; j++)
        !          1073:     {
        !          1074:       b = modii(mulii(a0, powmodulo(ginv0, n_q, p)), p);
        !          1075:       b = powmodulo(b, (GEN)qj[e-1-j], p);
        !          1076:       b = Fp_shanks(b, g_q, p, q);
        !          1077:       n_q = addii(n_q, mulii(b, (GEN)qj[j]));
        !          1078:     }
        !          1079:     v[i] = lmodulcp(n_q, (GEN)qj[e]);
        !          1080:   }
        !          1081:   return gerepileuptoint(av, lift(chinese(v,NULL)));
        !          1082: }
        !          1083:
        !          1084: /* smallest n >= 0 such that g0^n=x modulo pr, assume g0 reduced
        !          1085:  * q = order of g0  (Npr - 1 if q = NULL)
        !          1086:  * TODO: should be done in F_(p^f), not in Z_k mod pr (done for f=1) */
        !          1087: static GEN
        !          1088: nfshanks(GEN nf,GEN x,GEN g0,GEN pr,GEN prhall,GEN q)
        !          1089: {
        !          1090:   long av=avma,av1,lim,lbaby,i,k, f = itos((GEN)pr[4]);
        !          1091:   GEN p1,smalltable,giant,perm,v,g0inv,prh = (GEN)prhall[1];
        !          1092:   GEN multab, p = (GEN)pr[1];
        !          1093:
        !          1094:   x = lift_intern(nfreducemodpr(nf,x,prhall));
        !          1095:   p1 = q? q: addsi(-1, gpowgs(p,f));
        !          1096:   p1 = racine(p1);
        !          1097:   if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in nfshanks");
        !          1098:   lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
        !          1099:   g0inv = lift_intern(element_invmodpr(nf,g0,prhall));
        !          1100:   p1 = x;
        !          1101:
        !          1102:   multab = get_multab(nf, g0inv);
        !          1103:   for (i=lg(multab)-1; i; i--)
        !          1104:     multab[i]=(long)FpV_red((GEN)multab[i], p);
        !          1105:
        !          1106:   for (i=1;;i++)
        !          1107:   {
        !          1108:     if (isnfscalar(p1) && gcmp1((GEN)p1[1])) { avma=av; return stoi(i-1); }
        !          1109:
        !          1110:     smalltable[i]=(long)p1; if (i==lbaby) break;
        !          1111:     p1 = mul_matvec_mod_pr(multab,p1,prh);
        !          1112:   }
        !          1113:   giant=lift_intern(element_divmodpr(nf,x,p1,prhall));
        !          1114:   p1=cgetg(lbaby+1,t_VEC);
        !          1115:   perm = gen_sort(smalltable, cmp_IND | cmp_C, cmp_vecint);
        !          1116:   for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
        !          1117:   smalltable=p1; p1=giant;
        !          1118:
        !          1119:   multab = get_multab(nf, giant);
        !          1120:   for (i=lg(multab)-1; i; i--)
        !          1121:     multab[i]=(long)FpV_red((GEN)multab[i], p);
        !          1122:
        !          1123:   av1 = avma; lim=stack_lim(av1,2);
        !          1124:   for (k=1;;k++)
        !          1125:   {
        !          1126:     i=tablesearch(smalltable,p1,cmp_vecint);
        !          1127:     if (i)
        !          1128:     {
        !          1129:       v=addis(mulss(lbaby-1,k),perm[i]);
        !          1130:       return gerepileuptoint(av,addsi(-1,v));
        !          1131:     }
        !          1132:     p1 = mul_matvec_mod_pr(multab,p1,prh);
        !          1133:
        !          1134:     if (low_stack(lim, stack_lim(av1,2)))
        !          1135:     {
        !          1136:       if(DEBUGMEM>1) err(warnmem,"nfshanks");
        !          1137:       p1 = gerepileupto(av1, p1);
        !          1138:     }
        !          1139:   }
        !          1140: }
        !          1141:
        !          1142: /* same in nf.zk / pr
        !          1143:  * TODO: should be done in F_(p^f), not in Z_k mod pr (done for f=1) */
        !          1144: GEN
        !          1145: nf_PHlog(GEN nf, GEN a, GEN g, GEN pr, GEN prhall)
        !          1146: {
        !          1147:   ulong av = avma;
        !          1148:   GEN v,t0,a0,b,q,g_q,n_q,ginv0,qj,ginv,ord,fa,ex;
        !          1149:   long e,i,j,l;
        !          1150:
        !          1151:   a = lift_intern(nfreducemodpr(nf,a,prhall));
        !          1152:   if (isnfscalar(a))
        !          1153:   { /* can be done in Fp^* */
        !          1154:     GEN p = (GEN)pr[1], ordp = subis(p, 1);
        !          1155:     a = (GEN)a[1];
        !          1156:     if (gcmp1(a) || egalii(p, gdeux)) { avma = av; return gzero; }
        !          1157:     ord = subis(idealnorm(nf,pr), 1);
        !          1158:     if (egalii(a, ordp)) /* -1 */
        !          1159:       return gerepileuptoint(av, shifti(ord,-1));
        !          1160:
        !          1161:     if (egalii(ord, ordp))
        !          1162:       q = NULL;
        !          1163:     else /* we want < g > = Fp^* */
        !          1164:     {
        !          1165:       q = divii(ord,ordp);
        !          1166:       g = element_powmodpr(nf,g,q,prhall);
        !          1167:     }
        !          1168:     g = lift_intern((GEN)g[1]);
        !          1169:     n_q = Fp_PHlog(a,g,p,NULL);
        !          1170:     if (q) n_q = mulii(q, n_q);
        !          1171:     return gerepileuptoint(av, n_q);
        !          1172:   }
        !          1173:   ord = subis(idealnorm(nf,pr),1);
        !          1174:   fa = factor(ord); ex = (GEN)fa[2];
        !          1175:   fa = (GEN)fa[1];
        !          1176:   l = lg(fa);
        !          1177:   ginv = lift_intern(element_invmodpr(nf, g, prhall));
        !          1178:   v = cgetg(l, t_VEC);
        !          1179:   for (i=1; i<l; i++)
        !          1180:   {
        !          1181:     q = (GEN)fa[i];
        !          1182:     e = itos((GEN)ex[i]);
        !          1183:     if (DEBUGLEVEL>5)
        !          1184:       fprintferr("nf_Pohlig-Hellman: DL mod %Z^%ld\n",q,e);
        !          1185:     qj = new_chunk(e+1); qj[0] = un;
        !          1186:     for (j=1; j<=e; j++) qj[j] = lmulii((GEN)qj[j-1], q);
        !          1187:     t0 = divii(ord, (GEN)qj[e]);
        !          1188:     a0 = element_powmodpr(nf, a, t0, prhall);
        !          1189:     ginv0 = element_powmodpr(nf, ginv, t0, prhall); /* order q^e */
        !          1190:     g_q = element_powmodpr(nf, g, divii(ord,q), prhall); /* order q */
        !          1191:     n_q = gzero;
        !          1192:     for (j=0; j<e; j++)
        !          1193:     {
        !          1194:       b = element_mulmodpr(nf,a0,
        !          1195:                               element_powmodpr(nf, ginv0, n_q, prhall), prhall);
        !          1196:       b = element_powmodpr(nf, b, (GEN)qj[e-1-j], prhall);
        !          1197:       b = nfshanks(nf,b, g_q, pr,prhall, q);
        !          1198:       n_q = addii(n_q, mulii(b, (GEN)qj[j]));
        !          1199:     }
        !          1200:     v[i] = lmodulcp(n_q, (GEN)qj[e]);
        !          1201:   }
        !          1202:   return gerepileuptoint(av, lift(chinese(v,NULL)));
        !          1203: }
        !          1204:
        !          1205: GEN
        !          1206: znlog(GEN x, GEN g0)
        !          1207: {
        !          1208:   long av=avma;
        !          1209:   GEN p = (GEN)g0[1];
        !          1210:   if (typ(g0) != t_INTMOD)
        !          1211:     err(talker,"not an element of (Z/pZ)* in znlog");
        !          1212:   switch(typ(x))
        !          1213:   {
        !          1214:     case t_INT: break;
        !          1215:     default: x = gmul(x, gmodulsg(1,p));
        !          1216:     if (typ(x) != t_INTMOD)
        !          1217:       err(talker,"not an element of (Z/pZ)* in znlog");
        !          1218:     case t_INTMOD: x = (GEN)x[2]; break;
        !          1219:   }
        !          1220:   return gerepileuptoint(av, Fp_PHlog(x,(GEN)g0[2],p,NULL));
        !          1221: }
        !          1222:
        !          1223: GEN
        !          1224: dethnf(GEN mat)
        !          1225: {
        !          1226:   long i,l = lg(mat);
        !          1227:   ulong av;
        !          1228:   GEN s;
        !          1229:
        !          1230:   if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
        !          1231:   av = avma; s = gcoeff(mat,1,1);
        !          1232:   for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
        !          1233:   return av==avma? gcopy(s): gerepileupto(av,s);
        !          1234: }
        !          1235:
        !          1236: GEN
        !          1237: dethnf_i(GEN mat)
        !          1238: {
        !          1239:   long av,i,l = lg(mat);
        !          1240:   GEN s;
        !          1241:
        !          1242:   if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
        !          1243:   av = avma; s = gcoeff(mat,1,1);
        !          1244:   for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
        !          1245:   return gerepileuptoint(av,s);
        !          1246: }
        !          1247:
        !          1248: /* as above with cyc = diagonal(mat) */
        !          1249: GEN
        !          1250: detcyc(GEN cyc)
        !          1251: {
        !          1252:   long av,i,l = lg(cyc);
        !          1253:   GEN s;
        !          1254:
        !          1255:   if (l<3) return l<2? gun: icopy((GEN)cyc[1]);
        !          1256:   av = avma; s = (GEN)cyc[1];
        !          1257:   for (i=2; i<l; i++) s = mulii(s,(GEN)cyc[i]);
        !          1258:   return gerepileuptoint(av,s);
        !          1259: }
        !          1260:
        !          1261: /* (U,V) = 1. Return y = x mod U, = 1 mod V (uv: u + v = 1, u in U, v in V) */
        !          1262: static GEN
        !          1263: makeprimetoideal(GEN nf,GEN UV,GEN uv,GEN x)
        !          1264: {
        !          1265:   GEN y = gadd((GEN)uv[1], element_mul(nf,x,(GEN)uv[2]));
        !          1266:   return nfreducemodideal(nf,y,UV);
        !          1267: }
        !          1268:
        !          1269: static GEN
        !          1270: makeprimetoidealvec(GEN nf,GEN UV,GEN uv,GEN gen)
        !          1271: {
        !          1272:   long i, lx = lg(gen);
        !          1273:   GEN y = cgetg(lx,t_VEC);
        !          1274:   for (i=1; i<lx; i++)
        !          1275:     y[i] = (long)makeprimetoideal(nf,UV,uv,(GEN)gen[i]);
        !          1276:   return y;
        !          1277: }
        !          1278:
        !          1279: /* Given an ideal pr^ep, and an integral ideal x (in HNF form) compute a list
        !          1280:  * of vectors, each with 5 components as follows :
        !          1281:  * [[clh],[gen1],[gen2],[signat2],U.X^-1]. Each component corresponds to
        !          1282:  * d_i,g_i,g'_i,s_i.  Generators g_i are not necessarily prime to x, the
        !          1283:  * generators g'_i are. signat2 is the (horizontal) vector made of the
        !          1284:  * signatures (column vectors) of the g'_i. If x = NULL, the original ideal
        !          1285:  * was a prime power
        !          1286:  */
        !          1287: static GEN
        !          1288: zprimestar(GEN nf,GEN pr,GEN ep,GEN x,GEN arch)
        !          1289: {
        !          1290:   long av=avma,av1,N,f,nbp,j,n,m,tetpil,i,e,a,b;
        !          1291:   GEN prh,p,pf_1,list,v,p1,p3,p4,prk,uv,g0,newgen,pra,prb;
        !          1292:   GEN *gptr[2];
        !          1293:
        !          1294:   if(DEBUGLEVEL>3) { fprintferr("treating pr = %Z ^ %Z\n",pr,ep); flusherr(); }
        !          1295:   prh=prime_to_ideal(nf,pr); N=lg(prh)-1;
        !          1296:   f=itos((GEN)pr[4]); p=(GEN)pr[1];
        !          1297:
        !          1298:   pf_1 = addis(gpowgs(p,f), -1);
        !          1299:   v = zerocol(N);
        !          1300:   if (f==1) v[1]=gener(p)[2];
        !          1301:   else
        !          1302:   {
        !          1303:     GEN prhall = cgetg(3,t_VEC);
        !          1304:     long psim;
        !          1305:     if (is_bigint(p)) err(talker,"prime too big in zprimestar");
        !          1306:     psim = itos(p);
        !          1307:     list = (GEN)factor(pf_1)[1]; nbp=lg(list)-1;
        !          1308:     prhall[1]=(long)prh; prhall[2]=zero;
        !          1309:     for (n=psim; n>=0; n++)
        !          1310:     {
        !          1311:       m=n;
        !          1312:       for (i=1; i<=N; i++)
        !          1313:        if (!gcmp1(gcoeff(prh,i,i))) { v[i]=lstoi(m%psim); m/=psim; }
        !          1314:       for (j=1; j<=nbp; j++)
        !          1315:       {
        !          1316:         p1 = divii(pf_1,(GEN)list[j]);
        !          1317:        p1 = lift_intern(element_powmodpr(nf,v,p1,prhall));
        !          1318:         if (isnfscalar(p1) && gcmp1((GEN)p1[1])) break;
        !          1319:       }
        !          1320:       if (j>nbp) break;
        !          1321:     }
        !          1322:     if (n < 0) err(talker,"prime too big in zprimestar");
        !          1323:   }
        !          1324:   /* v generates  (Z_K / pr)^* */
        !          1325:   if(DEBUGLEVEL>3) {fprintferr("v computed\n");flusherr();}
        !          1326:   e = itos(ep); prk=(e==1)? pr: idealpow(nf,pr,ep);
        !          1327:   if(DEBUGLEVEL>3) {fprintferr("prk computed\n");flusherr();}
        !          1328:   g0 = v;
        !          1329:   uv = NULL; /* gcc -Wall */
        !          1330:   if (x)
        !          1331:   {
        !          1332:     uv = idealaddtoone(nf,prk, idealdivexact(nf,x,prk));
        !          1333:     g0 = makeprimetoideal(nf,x,uv,v);
        !          1334:     if(DEBUGLEVEL>3) {fprintferr("g0 computed\n");flusherr();}
        !          1335:   }
        !          1336:
        !          1337:   p1 = cgetg(6,t_VEC); list = _vec(p1);
        !          1338:   p1[1] = (long)_vec(pf_1);
        !          1339:   p1[2] = (long)_vec(v);
        !          1340:   p1[3] = (long)_vec(g0);
        !          1341:   p1[4] = (long)_vec(zsigne(nf,g0,arch));
        !          1342:   p1[5] = un;
        !          1343:   if (e==1) return gerepilecopy(av, list);
        !          1344:
        !          1345:   a=1; b=2; av1=avma;
        !          1346:   pra = prh; prb = (e==2)? prk: idealpow(nf,pr,gdeux);
        !          1347:   for(;;)
        !          1348:   {
        !          1349:     if(DEBUGLEVEL>3)
        !          1350:       {fprintferr("  treating a = %ld, b = %ld\n",a,b); flusherr();}
        !          1351:     p1 = zidealij(pra,prb);
        !          1352:     newgen = dummycopy((GEN)p1[2]);
        !          1353:     p3 = cgetg(lg(newgen),t_VEC);
        !          1354:     if(DEBUGLEVEL>3) {fprintferr("zidealij done\n"); flusherr();}
        !          1355:     for (i=1; i<lg(newgen); i++)
        !          1356:     {
        !          1357:       if (x) newgen[i]=(long)makeprimetoideal(nf,x,uv,(GEN)newgen[i]);
        !          1358:       p3[i]=(long)zsigne(nf,(GEN)newgen[i],arch);
        !          1359:     }
        !          1360:     p4 = cgetg(6,t_VEC);
        !          1361:     p4[1] = p1[1];
        !          1362:     p4[2] = p1[2];
        !          1363:     p4[3] = (long)newgen;
        !          1364:     p4[4] = (long)p3;
        !          1365:     p4[5] = p1[3]; p4 = _vec(p4);
        !          1366:
        !          1367:     a=b; b=min(e,b<<1); tetpil = avma;
        !          1368:     list = concat(list, p4);
        !          1369:     if (a==b) return gerepile(av,tetpil,list);
        !          1370:
        !          1371:     pra = gcopy(prb);
        !          1372:     gptr[0]=&pra; gptr[1]=&list;
        !          1373:     gerepilemanysp(av1,tetpil,gptr,2);
        !          1374:     prb = (b==(a<<1))? idealpow(nf,pra,gdeux): prk;
        !          1375:   }
        !          1376: }
        !          1377:
        !          1378: /* x ideal, compute elements in 1+x whose sign matrix is invertible */
        !          1379: GEN
        !          1380: zarchstar(GEN nf,GEN x,GEN arch,long nba)
        !          1381: {
        !          1382:   long av,N,i,j,k,r,rr,limr,zk,lgmat;
        !          1383:   GEN p1,y,bas,genarch,mat,lambda,nfun,vun;
        !          1384:
        !          1385:   y = cgetg(4,t_VEC);
        !          1386:   if (!nba)
        !          1387:   {
        !          1388:     y[1]=lgetg(1,t_VEC);
        !          1389:     y[2]=lgetg(1,t_VEC);
        !          1390:     y[3]=lgetg(1,t_MAT); return y;
        !          1391:   }
        !          1392:   N = degpol(nf[1]);
        !          1393:   p1 = cgetg(nba+1,t_VEC); for (i=1; i<=nba; i++) p1[i] = deux;
        !          1394:   y[1] = (long)p1; av = avma;
        !          1395:   if (N == 1)
        !          1396:   {
        !          1397:     p1 = gerepileuptoint(av, subsi(1, shifti(gcoeff(x,1,1),1)));
        !          1398:     y[2] = (long)_vec(_col(p1));
        !          1399:     y[3] = (long)gscalmat(gun,1); return y;
        !          1400:   }
        !          1401:   zk = ideal_is_zk(x,N);
        !          1402:   genarch = cgetg(nba+1,t_VEC);
        !          1403:   mat = cgetg(nba+1,t_MAT); setlg(mat,2);
        !          1404:   for (i=1; i<=nba; i++) mat[i] = lgetg(nba+1, t_MAT);
        !          1405:   lambda = cgetg(N+1, t_VECSMALL);
        !          1406:
        !          1407:   bas = dummycopy(gmael(nf,5,1)); r = lg(arch);
        !          1408:   for (k=1,i=1; i<r; i++)
        !          1409:     if (signe(arch[i]))
        !          1410:     {
        !          1411:       if (k < i)
        !          1412:         for (j=1; j<=N; j++) coeff(bas,k,j) = coeff(bas,i,j);
        !          1413:       k++;
        !          1414:     }
        !          1415:   r = nba+1; for (i=1; i<=N; i++) setlg(bas[i], r);
        !          1416:   if (!zk)
        !          1417:   {
        !          1418:     x = gmul(x,lllintpartial(x));
        !          1419:     bas = gmul(bas, x);
        !          1420:   }
        !          1421:
        !          1422:   vun = cgetg(nba+1, t_COL);
        !          1423:   for (i=1; i<=nba; i++) vun[i] = un;
        !          1424:   nfun = gscalcol(gun, N); lgmat = 1;
        !          1425:   for (r=1, rr=3; ; r<<=1, rr=(r<<1)+1)
        !          1426:   {
        !          1427:     if (DEBUGLEVEL>5) fprintferr("zarchstar: r = %ld\n",r);
        !          1428:     p1 = gpowgs(stoi(rr),N);
        !          1429:     limr = is_bigint(p1)? BIGINT: p1[2];
        !          1430:     limr = (limr-1)>>1;
        !          1431:     k = zk? rr: 1; /* to avoid 0 */
        !          1432:     for ( ; k<=limr; k++)
        !          1433:     {
        !          1434:       long av1=avma, kk = k;
        !          1435:       GEN alpha = vun;
        !          1436:       for (i=1; i<=N; i++)
        !          1437:       {
        !          1438:         lambda[i] = (kk+r)%rr - r; kk/=rr;
        !          1439:         if (lambda[i])
        !          1440:           alpha = gadd(alpha, gmulsg(lambda[i],(GEN)bas[i]));
        !          1441:       }
        !          1442:       p1 = (GEN)mat[lgmat];
        !          1443:       for (i=1; i<=nba; i++)
        !          1444:         p1[i] = (gsigne((GEN)alpha[i]) > 0)? zero: un;
        !          1445:
        !          1446:       if (ker_trivial_mod_p(mat, gdeux)) avma = av1;
        !          1447:       else
        !          1448:       { /* new vector indep. of previous ones */
        !          1449:         avma = av1; alpha = nfun;
        !          1450:         for (i=1; i<=N; i++)
        !          1451:           if (lambda[i])
        !          1452:             alpha = gadd(alpha, gmulsg(lambda[i],(GEN)x[i]));
        !          1453:        genarch[lgmat++] = (long)alpha;
        !          1454:        if (lgmat > nba)
        !          1455:        {
        !          1456:           GEN *gptr[2];
        !          1457:           mat = ginv(mat); gptr[0]=&genarch; gptr[1]=&mat;
        !          1458:           gerepilemany(av,gptr,2);
        !          1459:          y[2]=(long)genarch;
        !          1460:           y[3]=(long)mat; return y;
        !          1461:        }
        !          1462:         setlg(mat,lgmat+1);
        !          1463:       }
        !          1464:     }
        !          1465:   }
        !          1466: }
        !          1467:
        !          1468: GEN
        !          1469: zinternallog_pk(GEN nf, GEN a0, GEN y, GEN pr, GEN prk, GEN list, GEN *psigne)
        !          1470: {
        !          1471:   GEN a = a0, L,e,p1,cyc,gen;
        !          1472:   long i,j, llist = lg(list)-1;
        !          1473:
        !          1474:   e = gzero; /* gcc -Wall */
        !          1475:   for (j=1; j<=llist; j++)
        !          1476:   {
        !          1477:     L = (GEN)list[j]; cyc=(GEN)L[1]; gen=(GEN)L[2];
        !          1478:     if (j==1)
        !          1479:       p1 = nf_PHlog(nf,a,(GEN)gen[1],pr, nfmodprinit(nf,pr));
        !          1480:     else
        !          1481:     {
        !          1482:       p1 = (GEN)a[1]; a[1] = laddsi(-1,(GEN)a[1]);
        !          1483:       e = gmul((GEN)L[5],a); a[1] = (long)p1;
        !          1484:       /* here lg(e) == lg(cyc) */
        !          1485:       p1 = (GEN)e[1];
        !          1486:     }
        !          1487:     for (i=1; i<lg(cyc); i++, p1 = (GEN)e[i])
        !          1488:     {
        !          1489:       p1 = modii(negi(p1), (GEN)cyc[i]);
        !          1490:       *++y = lnegi(p1);
        !          1491:       if (!signe(p1)) continue;
        !          1492:
        !          1493:       if (psigne && mpodd(p1)) *psigne = gadd(*psigne,gmael(L,4,i));
        !          1494:       if (j == llist) continue;
        !          1495:
        !          1496:       if (DEBUGLEVEL>5) fprintferr("do element_powmodideal\n");
        !          1497:       p1 = element_powmodideal(nf,(GEN)gen[i],p1,prk);
        !          1498:       a = element_mulmodideal(nf,a,p1,prk);
        !          1499:     }
        !          1500:   }
        !          1501:   return y;
        !          1502: }
        !          1503:
        !          1504: /* Retourne la decomposition de a sur les nbgen generateurs successifs
        !          1505:  * contenus dans list_set et si index !=0 on ne fait ce calcul que pour
        !          1506:  * l'ideal premier correspondant a cet index en mettant a 0 les autres
        !          1507:  * composantes
        !          1508:  */
        !          1509: static GEN
        !          1510: zinternallog(GEN nf,GEN a,GEN list_set,long nbgen,GEN arch,GEN fa,long index)
        !          1511: {
        !          1512:   GEN prlist,ep,y0,y,ainit,list,pr,prk,cyc,psigne,p1,p2;
        !          1513:   long av,nbp,i,j,k;
        !          1514:
        !          1515:   y0 = y = cgetg(nbgen+1,t_COL); av=avma;
        !          1516:   prlist=(GEN)fa[1]; ep=(GEN)fa[2]; nbp=lg(ep)-1;
        !          1517:   i=typ(a); if (is_extscalar_t(i)) a = algtobasis(nf,a);
        !          1518:   if (DEBUGLEVEL>3)
        !          1519:   {
        !          1520:     fprintferr("entering zinternallog, "); flusherr();
        !          1521:     if (DEBUGLEVEL>5) fprintferr("with a = %Z\n",a);
        !          1522:   }
        !          1523:   ainit = a; psigne = zsigne(nf,ainit,arch);
        !          1524:   p2 = NULL; /* gcc -Wall */
        !          1525:   for (k=1; k<=nbp; k++)
        !          1526:   {
        !          1527:     list = (GEN)list_set[k];
        !          1528:     if (index && index!=k)
        !          1529:     {
        !          1530:       for (j=1; j<lg(list); j++)
        !          1531:       {
        !          1532:         cyc = gmael(list,j,1);
        !          1533:         for (i=1; i<lg(cyc); i++) *++y = zero;
        !          1534:       }
        !          1535:       continue;
        !          1536:     }
        !          1537:     pr = (GEN)prlist[k]; prk = idealpow(nf,pr,(GEN)ep[k]);
        !          1538:     y = zinternallog_pk(nf, ainit, y, pr, prk, list, &psigne);
        !          1539:   }
        !          1540:   p1 = lift_intern(gmul(gmael(list_set,k,3), psigne));
        !          1541:   avma=av; for (i=1; i<lg(p1); i++) *++y = p1[i];
        !          1542:   if (DEBUGLEVEL>3) { fprintferr("leaving\n"); flusherr(); }
        !          1543:   for (i=1; i<=nbgen; i++) y0[i] = licopy((GEN)y0[i]);
        !          1544:   return y0;
        !          1545: }
        !          1546:
        !          1547: static GEN
        !          1548: compute_gen(GEN nf, GEN u1, GEN gen, GEN bid)
        !          1549: {
        !          1550:   long i, c = lg(u1);
        !          1551:   GEN basecl = cgetg(c,t_VEC);
        !          1552:
        !          1553:   for (i=1; i<c; i++)
        !          1554:     basecl[i] = (long)famat_to_nf_modidele(nf, gen, (GEN)u1[i], bid);
        !          1555:   return basecl;
        !          1556: }
        !          1557:
        !          1558: /* Compute [[ideal,arch], [h,[cyc],[gen]], idealfact, [liste], U]
        !          1559:    gen not computed unless add_gen != 0 */
        !          1560: GEN
        !          1561: zidealstarinitall(GEN nf, GEN ideal,long add_gen)
        !          1562: {
        !          1563:   long av = avma,i,j,k,nba,nbp,R1,nbgen,cp;
        !          1564:   GEN p1,y,h,cyc,U,u1,fa,fa2,ep,x,arch,list,sarch,gen;
        !          1565:
        !          1566:   nf = checknf(nf); R1=itos(gmael(nf,2,1));
        !          1567:   if (typ(ideal)==t_VEC && lg(ideal)==3)
        !          1568:   {
        !          1569:     arch=(GEN)ideal[2]; ideal = (GEN)ideal[1];
        !          1570:     i = typ(arch);
        !          1571:     if (!is_vec_t(i) || lg(arch) != R1+1)
        !          1572:       err(talker,"incorrect archimedean component in zidealstarinit");
        !          1573:     for (nba=0,i=1; i<=R1; i++)
        !          1574:       if (signe(arch[i])) nba++;
        !          1575:   }
        !          1576:   else
        !          1577:   {
        !          1578:     arch=cgetg(R1+1,t_VEC);
        !          1579:     for (nba=0,i=1; i<=R1; i++) arch[i]=zero;
        !          1580:   }
        !          1581:   x = idealhermite(nf,ideal);
        !          1582:   if (!gcmp1(denom(x)))
        !          1583:     err(talker,"zidealstarinit needs an integral ideal: %Z",x);
        !          1584:   p1=cgetg(3,t_VEC); ideal=p1;
        !          1585:   p1[1]=(long)x;
        !          1586:   p1[2]=(long)arch;
        !          1587:
        !          1588:   fa=idealfactor(nf,x); list=(GEN)fa[1]; ep=(GEN)fa[2];
        !          1589:   nbp=lg(list)-1; fa2=cgetg(nbp+2,t_VEC);
        !          1590:
        !          1591:   gen = cgetg(1,t_VEC);
        !          1592:   p1 = (nbp==1)? (GEN)NULL: x;
        !          1593:   for (i=1; i<=nbp; i++)
        !          1594:   {
        !          1595:     GEN L = zprimestar(nf,(GEN)list[i],(GEN)ep[i],p1,arch);
        !          1596:     fa2[i] = (long)L;
        !          1597:     for (j=1; j<lg(L); j++) gen = concatsp(gen,gmael(L,j,3));
        !          1598:   }
        !          1599:   sarch = zarchstar(nf,x,arch,nba);
        !          1600:   fa2[i]=(long)sarch;
        !          1601:   gen = concatsp(gen,(GEN)sarch[2]);
        !          1602:
        !          1603:   nbgen = lg(gen)-1;
        !          1604:   h=cgetg(nbgen+1,t_MAT); cp=0;
        !          1605:   for (i=1; i<=nbp; i++)
        !          1606:   {
        !          1607:     list = (GEN)fa2[i];
        !          1608:     for (j=1; j<lg(list); j++)
        !          1609:     {
        !          1610:       GEN a, L = (GEN)list[j], e = (GEN)L[1], g = (GEN)L[3];
        !          1611:       for (k=1; k<lg(g); k++)
        !          1612:       {
        !          1613:        a = element_powmodidele(nf,(GEN)g[k],(GEN)e[k],ideal,sarch);
        !          1614:        h[++cp] = lneg(zinternallog(nf,a,fa2,nbgen,arch,fa,i));
        !          1615:        coeff(h,cp,cp) = e[k];
        !          1616:       }
        !          1617:     }
        !          1618:   }
        !          1619:   for (j=1; j<=nba; j++) { h[++cp]=(long)zerocol(nbgen); coeff(h,cp,cp)=deux; }
        !          1620:   if (cp!=nbgen) err(talker,"bug in zidealstarinit");
        !          1621:   h = hnfall_i(h,NULL,0);
        !          1622:   cyc = smithrel(h, &U, add_gen? &u1: NULL);
        !          1623:   p1 = cgetg(add_gen? 4: 3, t_VEC);
        !          1624:   p1[1] = (long)detcyc(cyc);
        !          1625:   p1[2] = (long)cyc;
        !          1626:
        !          1627:   y = cgetg(6,t_VEC);
        !          1628:   y[1] = (long)ideal;
        !          1629:   y[2] = (long)p1;
        !          1630:   y[3] = (long)fa;
        !          1631:   y[4] = (long)fa2;
        !          1632:   y[5] = (long)U;
        !          1633:   if (add_gen) p1[3] = (long)compute_gen(nf,u1,gen, y);
        !          1634:   return gerepilecopy(av, y);
        !          1635: }
        !          1636:
        !          1637: GEN
        !          1638: zidealstarinitgen(GEN nf, GEN ideal)
        !          1639: {
        !          1640:   return zidealstarinitall(nf,ideal,1);
        !          1641: }
        !          1642:
        !          1643: GEN
        !          1644: zidealstarinit(GEN nf, GEN ideal)
        !          1645: {
        !          1646:   return zidealstarinitall(nf,ideal,0);
        !          1647: }
        !          1648:
        !          1649: GEN
        !          1650: zidealstar(GEN nf, GEN ideal)
        !          1651: {
        !          1652:   long av = avma;
        !          1653:   GEN y = zidealstarinitall(nf,ideal,1);
        !          1654:   return gerepilecopy(av,(GEN)y[2]);
        !          1655: }
        !          1656:
        !          1657: GEN
        !          1658: idealstar0(GEN nf, GEN ideal,long flag)
        !          1659: {
        !          1660:   switch(flag)
        !          1661:   {
        !          1662:     case 0: return zidealstar(nf,ideal);
        !          1663:     case 1: return zidealstarinit(nf,ideal);
        !          1664:     case 2: return zidealstarinitgen(nf,ideal);
        !          1665:     default: err(flagerr,"idealstar");
        !          1666:   }
        !          1667:   return NULL; /* not reached */
        !          1668: }
        !          1669:
        !          1670: void
        !          1671: check_nfelt(GEN x, GEN *den)
        !          1672: {
        !          1673:   long l = lg(x), i;
        !          1674:   GEN t, d = NULL;
        !          1675:   if (typ(x) != t_COL) err(talker,"%Z not a nfelt", x);
        !          1676:   for (i=1; i<l; i++)
        !          1677:   {
        !          1678:     t = (GEN)x[i];
        !          1679:     switch (typ(t))
        !          1680:     {
        !          1681:       case t_INT: case t_INTMOD: break;
        !          1682:       case t_FRACN: case t_FRAC:
        !          1683:         if (!d) d = (GEN)t[2]; else d = mpppcm(d, (GEN)t[2]);
        !          1684:         break;
        !          1685:       default: err(talker,"%Z not a nfelt", x);
        !          1686:     }
        !          1687:   }
        !          1688:   *den = d;
        !          1689: }
        !          1690:
        !          1691: /* Given x (not necessarily integral), and bid as output by zidealstarinit,
        !          1692:  * compute the vector of components on the generators bid[2].
        !          1693:  * Assume (x,bid) = 1 */
        !          1694: GEN
        !          1695: zideallog(GEN nf, GEN x, GEN bid)
        !          1696: {
        !          1697:   long av,l,i,N,c;
        !          1698:   GEN fa,fa2,ideal,arch,den,p1,cyc,y;
        !          1699:
        !          1700:   nf=checknf(nf); checkbid(bid);
        !          1701:   cyc=gmael(bid,2,2); c=lg(cyc);
        !          1702:   y=cgetg(c,t_COL); av=avma;
        !          1703:   N = degpol(nf[1]); ideal = (GEN) bid[1];
        !          1704:   if (typ(ideal)==t_VEC && lg(ideal)==3)
        !          1705:     arch = (GEN)ideal[2];
        !          1706:   else
        !          1707:     arch = NULL;
        !          1708:   switch(typ(x))
        !          1709:   {
        !          1710:     case t_INT: case t_FRAC: case t_FRACN:
        !          1711:       x = gscalcol_i(x,N); break;
        !          1712:     case t_POLMOD: case t_POL:
        !          1713:       x = algtobasis(nf,x); break;
        !          1714:     case t_MAT:
        !          1715:       return famat_ideallog(nf,(GEN)x[1],(GEN)x[2],bid);
        !          1716:     case t_COL: break;
        !          1717:     default: err(talker,"not an element in zideallog");
        !          1718:   }
        !          1719:   if (lg(x) != N+1) err(talker,"not an element in zideallog");
        !          1720:   check_nfelt(x, &den);
        !          1721:   if (den)
        !          1722:   {
        !          1723:     GEN g = cgetg(3, t_COL);
        !          1724:     GEN e = cgetg(3, t_COL);
        !          1725:     g[1] = lmul(x,den); e[1] = un;
        !          1726:     g[2] = (long)den;   e[2] = lstoi(-1);
        !          1727:     p1 = famat_ideallog(nf, g, e, bid);
        !          1728:   }
        !          1729:   else
        !          1730:   {
        !          1731:     l=lg(bid[5])-1; fa=(GEN)bid[3]; fa2=(GEN)bid[4];
        !          1732:     p1 = zinternallog(nf,x,fa2,l,arch,fa,0);
        !          1733:     p1 = gmul((GEN)bid[5],p1); /* apply smith */
        !          1734:   }
        !          1735:   if (lg(p1)!=c) err(bugparier,"zideallog");
        !          1736:   for (i=1; i<c; i++)
        !          1737:     y[i] = lmodii((GEN)p1[i],(GEN)cyc[i]);
        !          1738:   avma=av; /* following line does a gerepile ! */
        !          1739:   for (i=1; i<c; i++)
        !          1740:     y[i] = (long)icopy((GEN)y[i]);
        !          1741:   return y;
        !          1742: }
        !          1743:
        !          1744: /* bid1, bid2: output from 'zidealstarinit' for coprime modules m1 and m2
        !          1745:  * (without arch. part).
        !          1746:  * Output: zidealstarinit [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U]
        !          1747:  * associated to m1 m2 */
        !          1748: GEN
        !          1749: zidealstarinitjoin(GEN nf, GEN bid1, GEN bid2, long add_gen)
        !          1750: {
        !          1751:   long av=avma,i,nbp,nbgen,lx1,lx2,l1,l2,lx;
        !          1752:   GEN module1,module2,struct1,struct2,fact1,fact2,liste1,liste2;
        !          1753:   GEN module,liste,fact,U,U1,U2,cyc,cyc1,cyc2,uv;
        !          1754:   GEN p1,p2,y,u1,x,fa1,fa2, gen = add_gen? gun: NULL;
        !          1755:
        !          1756:   nf = checknf(nf); checkbid(bid1); checkbid(bid2);
        !          1757:   module1 = (GEN)bid1[1]; struct1 = (GEN)bid1[2]; fact1 = (GEN)bid1[3];
        !          1758:   module2 = (GEN)bid2[1]; struct2 = (GEN)bid2[2]; fact2 = (GEN)bid2[3];
        !          1759:   x = idealmul(nf,(GEN)module1[1],(GEN)module2[1]);
        !          1760:   module = cgetg(3,t_VEC);
        !          1761:   module[1] = (long)x;
        !          1762:   module[2] = module1[2];
        !          1763:   if (!gcmp0((GEN)module1[2]) || !gcmp0((GEN)module2[2]))
        !          1764:     err(talker,"non-0 Archimedian components in zidealstarinitjoin");
        !          1765:
        !          1766:   fa1 = (GEN)fact1[1]; lx1 = lg(fa1);
        !          1767:   fa2 = (GEN)fact2[1]; lx2 = lg(fa2);
        !          1768:   fact = cgetg(3,t_MAT);
        !          1769:   fact[1] = (long)concatsp(fa1,fa2);
        !          1770:   fact[2] = (long)concatsp((GEN)fact1[2],(GEN)fact2[2]);
        !          1771:   for (i=1; i<lx1; i++)
        !          1772:     if (isinvector(fa2,(GEN)fa1[i],lx2-1))
        !          1773:       err(talker,"noncoprime ideals in zidealstarinitjoin");
        !          1774:
        !          1775:   liste1 = (GEN)bid1[4]; lx1 = lg(liste1);
        !          1776:   liste2 = (GEN)bid2[4]; lx2 = lg(liste2);
        !          1777:   /* concat (liste1 - last elt [zarchstar]) + liste2 */
        !          1778:   lx = lx1+lx2-2; liste = cgetg(lx,t_VEC);
        !          1779:   for (i=1; i<lx1-1; i++) liste[i] = liste1[i];
        !          1780:   for (   ; i<lx; i++)    liste[i] = liste2[i-lx1+2];
        !          1781:
        !          1782:   U1 = (GEN)bid1[5]; cyc1 = (GEN)struct1[2]; l1 = lg(cyc1);
        !          1783:   U2 = (GEN)bid2[5]; cyc2 = (GEN)struct2[2]; l2 = lg(cyc2);
        !          1784:   nbgen = l1+l2-2;
        !          1785:   cyc = diagonal(concatsp(cyc1,cyc2));
        !          1786:   cyc = smithrel(cyc, &U, gen? &u1: NULL);
        !          1787:
        !          1788:   if (nbgen)
        !          1789:     U = concatsp(
        !          1790:       gmul(vecextract_i(U, 1,   l1-1), U1) ,
        !          1791:       gmul(vecextract_i(U, l1, nbgen), U2)
        !          1792:     );
        !          1793:
        !          1794:   if (gen)
        !          1795:   {
        !          1796:     if (lg(struct1)<=3 || lg(struct2)<=3)
        !          1797:       err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
        !          1798:
        !          1799:     uv = idealaddtoone(nf,(GEN)module1[1],(GEN)module2[1]);
        !          1800:     p1 = makeprimetoidealvec(nf,x,uv,(GEN)struct1[3]);
        !          1801:     p2 = (GEN)uv[1]; uv[1] = uv[2]; uv[2] = (long)p2;
        !          1802:     p2 = makeprimetoidealvec(nf,x,uv,(GEN)struct2[3]);
        !          1803:     gen = concatsp(p1,p2);
        !          1804:     nbp = lg((GEN)fact[1])-1;
        !          1805:   }
        !          1806:   p1 = cgetg(gen? 4: 3,t_VEC);
        !          1807:   p1[1] = (long)detcyc(cyc);
        !          1808:   p1[2] = (long)cyc;
        !          1809:
        !          1810:   y = cgetg(6,t_VEC);
        !          1811:   y[1] = (long)module;
        !          1812:   y[2] = (long)p1;
        !          1813:   y[3] = (long)fact;
        !          1814:   y[4] = (long)liste;
        !          1815:   y[5] = (long)U;
        !          1816:   if (gen) p1[3] = (long)compute_gen(nf,u1,gen, y);
        !          1817:   return gerepilecopy(av,y);
        !          1818: }
        !          1819:
        !          1820: /* bid1: output from 'zidealstarinit' for module m1 (without arch. part)
        !          1821:  * arch: archimedean part.
        !          1822:  * Output: zidealstarinit [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U]
        !          1823:  * associated to m1.arch */
        !          1824: GEN
        !          1825: zidealstarinitjoinarch(GEN nf, GEN bid1, GEN arch, long nba, long add_gen)
        !          1826: {
        !          1827:   long av=avma,i,nbp,lx1;
        !          1828:   GEN module1,struct1,fact1,liste1,U1,U;
        !          1829:   GEN module,liste,cyc,p1,y,u1,x,sarch, gen = add_gen? gun: NULL;
        !          1830:
        !          1831:   nf = checknf(nf); checkbid(bid1);
        !          1832:   module1 = (GEN)bid1[1]; struct1 = (GEN)bid1[2]; fact1 = (GEN)bid1[3];
        !          1833:   x = (GEN)module1[1];
        !          1834:   module = cgetg(3,t_VEC);
        !          1835:   module[1] = (long)x;
        !          1836:   module[2] = (long)arch;
        !          1837:   if (!gcmp0((GEN)module1[2]))
        !          1838:     err(talker,"non-0 Archimedian components in zidealstarinitjoinarch");
        !          1839:
        !          1840:   sarch = zarchstar(nf,x,arch,nba);
        !          1841:   liste1 = (GEN)bid1[4]; lx1 = lg(liste1);
        !          1842:   liste = cgetg(lx1,t_VEC);
        !          1843:   for (i=1; i<lx1-1; i++) liste[i]=liste1[i];
        !          1844:   liste[i] = (long)sarch;
        !          1845:
        !          1846:   U1 = (GEN)bid1[5]; lx1 = lg(U1);
        !          1847:   cyc = diagonal(concatsp((GEN)struct1[2],(GEN)sarch[1]));
        !          1848:   cyc = smithrel(cyc, &U, gen? &u1: NULL);
        !          1849:
        !          1850:   if (gen)
        !          1851:   {
        !          1852:     if (lg(struct1)<=3)
        !          1853:       err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
        !          1854:     gen = concatsp((GEN)struct1[3],(GEN)sarch[2]);
        !          1855:     nbp = lg((GEN)fact1[1])-1;
        !          1856:   }
        !          1857:   p1 = cgetg(gen? 4: 3, t_VEC);
        !          1858:   p1[1] = (long)detcyc(cyc);
        !          1859:   p1[2] = (long)cyc;
        !          1860:
        !          1861:   y = cgetg(6,t_VEC);
        !          1862:   y[1] = (long)module;
        !          1863:   y[2] = (long)p1;
        !          1864:   y[3] = (long)fact1;
        !          1865:   y[4] = (long)liste;
        !          1866:   y[5] = (long)U;
        !          1867:   if (gen) p1[3] = (long)compute_gen(nf,u1,gen, y);
        !          1868:   return gerepilecopy(av,y);
        !          1869: }
        !          1870:
        !          1871: /* calcule la matrice des zinternallog des unites */
        !          1872: GEN
        !          1873: logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid)
        !          1874: {
        !          1875:   long R,j,sizeh;
        !          1876:   GEN m,fa2,fa,arch;
        !          1877:
        !          1878:   R=lg(funits)-1; m=cgetg(R+2,t_MAT);
        !          1879:   fa2=(GEN)bid[4]; sizeh=lg(bid[5])-1; arch=gmael(bid,1,2);
        !          1880:   fa=(GEN)bid[3];
        !          1881:   m[1]=(long)zinternallog(nf,racunit,fa2,sizeh,arch,fa,0);
        !          1882:   for (j=2; j<=R+1; j++)
        !          1883:     m[j]=(long)zinternallog(nf,(GEN)funits[j-1],fa2,sizeh,arch,fa,0);
        !          1884:   return m;
        !          1885: }
        !          1886:
        !          1887: /* calcule la matrice des zinternallog des unites */
        !          1888: static GEN
        !          1889: logunitmatrixarch(GEN nf,GEN funits,GEN racunit,GEN bid)
        !          1890: {
        !          1891:   long R,j;
        !          1892:   GEN m,liste,structarch,arch;
        !          1893:
        !          1894:   R=lg(funits)-1; m=cgetg(R+2,t_MAT); arch=gmael(bid,1,2);
        !          1895:   liste=(GEN)bid[4]; structarch=(GEN)liste[lg(liste)-1];
        !          1896:   m[1]=(long)zsigne(nf,racunit,arch);
        !          1897:   for (j=2; j<=R+1; j++)
        !          1898:     m[j]=(long)zsigne(nf,(GEN)funits[j-1],arch);
        !          1899:   return lift_intern(gmul((GEN)structarch[3],m));
        !          1900: }
        !          1901:
        !          1902: static void
        !          1903: init_units(GEN bnf, GEN *funits, GEN *racunit)
        !          1904: {
        !          1905:   GEN p1;
        !          1906:   bnf = checkbnf(bnf); p1=(GEN)bnf[8];
        !          1907:   if (lg(p1)==5) *funits=(GEN)buchfu(bnf)[1];
        !          1908:   else
        !          1909:   {
        !          1910:     if (lg(p1)!=7) err(talker,"incorrect big number field");
        !          1911:     *funits=(GEN)p1[5];
        !          1912:   }
        !          1913:   *racunit=gmael(p1,4,2);
        !          1914: }
        !          1915:
        !          1916: /*  flag &1 : generateurs, sinon non
        !          1917:  *  flag &2 : unites, sinon pas.
        !          1918:  *  flag &4 : ideaux, sinon zidealstar.
        !          1919:  */
        !          1920: static GEN
        !          1921: ideallistzstarall(GEN bnf,long bound,long flag)
        !          1922: {
        !          1923:   byteptr ptdif=diffptr;
        !          1924:   long lim,av0=avma,av,i,j,k,l,q2,lp1,q;
        !          1925:   long do_gen = flag & 1, do_units = flag & 2, big_id = !(flag & 4);
        !          1926:   GEN y,nf,p,z,z2,p1,p2,p3,fa,pr,ideal,lu,lu2,funits,racunit,embunit;
        !          1927:
        !          1928:   nf = checknf(bnf);
        !          1929:   if (bound <= 0) return cgetg(1,t_VEC);
        !          1930:   z = cgetg(bound+1,t_VEC);
        !          1931:   for (i=2; i<=bound; i++) z[i] = lgetg(1,t_VEC);
        !          1932:
        !          1933:   ideal = idmat(degpol(nf[1]));
        !          1934:   if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
        !          1935:   z[1] = (long)_vec(ideal);
        !          1936:   if (do_units)
        !          1937:   {
        !          1938:     init_units(bnf,&funits,&racunit);
        !          1939:     lu = cgetg(bound+1,t_VEC);
        !          1940:     for (i=2; i<=bound; i++) lu[i]=lgetg(1,t_VEC);
        !          1941:     lu[1] = (long)_vec(logunitmatrix(nf,funits,racunit,ideal));
        !          1942:   }
        !          1943:
        !          1944:   p=cgeti(3); p[1]=evalsigne(1) | evallgefint(3);
        !          1945:   av=avma; lim=stack_lim(av,1);
        !          1946:   lu2 = embunit = NULL; /* gcc -Wall */
        !          1947:   if (bound > (long)maxprime()) err(primer1);
        !          1948:   for (p[2]=0; p[2]<=bound; )
        !          1949:   {
        !          1950:     p[2] += *ptdif++;
        !          1951:     if (DEBUGLEVEL>1) { fprintferr("%ld ",p[2]); flusherr(); }
        !          1952:     fa = primedec(nf,p);
        !          1953:     for (j=1; j<lg(fa); j++)
        !          1954:     {
        !          1955:       pr = (GEN)fa[j]; p1 = powgi(p,(GEN)pr[4]);
        !          1956:       if (is_bigint(p1) || (q = itos(p1)) > bound) continue;
        !          1957:
        !          1958:       q2=q; ideal=pr; z2=dummycopy(z);
        !          1959:       if (do_units) lu2=dummycopy(lu);
        !          1960:       for (l=2; ;l++)
        !          1961:       {
        !          1962:         if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
        !          1963:         if (do_units) embunit = logunitmatrix(nf,funits,racunit,ideal);
        !          1964:         for (i=q; i<=bound; i+=q)
        !          1965:         {
        !          1966:           p1 = (GEN)z[i/q]; lp1 = lg(p1);
        !          1967:           if (lp1 == 1) continue;
        !          1968:
        !          1969:           p2 = cgetg(lp1,t_VEC);
        !          1970:           for (k=1; k<lp1; k++)
        !          1971:             if (big_id)
        !          1972:               p2[k] = (long)zidealstarinitjoin(nf,(GEN)p1[k],ideal,do_gen);
        !          1973:             else
        !          1974:               p2[k] = (long)idealmul(nf,(GEN)p1[k],ideal);
        !          1975:           z2[i] = (long)concatsp((GEN)z2[i],p2);
        !          1976:           if (do_units)
        !          1977:           {
        !          1978:             p1 = (GEN)lu[i/q];
        !          1979:             p2 = cgetg(lp1,t_VEC);
        !          1980:             for (k=1; k<lp1; k++)
        !          1981:               p2[k] = (long)vconcat((GEN)p1[k],embunit);
        !          1982:             lu2[i] = (long)concatsp((GEN)lu2[i],p2);
        !          1983:           }
        !          1984:         }
        !          1985:         q *= q2; if ((ulong)q > (ulong)bound) break;
        !          1986:         ideal = idealpows(nf,pr,l);
        !          1987:       }
        !          1988:       z = z2; if (do_units) lu = lu2;
        !          1989:     }
        !          1990:     if (low_stack(lim, stack_lim(av,1)))
        !          1991:     {
        !          1992:       GEN *gptr[2]; gptr[0]=&z; gptr[1]=&lu;
        !          1993:       if(DEBUGMEM>1) err(warnmem,"ideallistzstarall");
        !          1994:       gerepilemany(av,gptr,do_units?2:1);
        !          1995:     }
        !          1996:   }
        !          1997:   if (!do_units) return gerepilecopy(av0, z);
        !          1998:   y = cgetg(3,t_VEC);
        !          1999:   y[1] = lcopy(z);
        !          2000:   lu2 = cgetg(lg(z),t_VEC);
        !          2001:   for (i=1; i<lg(z); i++)
        !          2002:   {
        !          2003:     p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
        !          2004:     p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
        !          2005:     for (j=1; j<lp1; j++) p3[j] = lmul(gmael(p1,j,5),(GEN)p2[j]);
        !          2006:   }
        !          2007:   y[2]=(long)lu2; return gerepileupto(av0, y);
        !          2008: }
        !          2009:
        !          2010: GEN
        !          2011: ideallist0(GEN bnf,long bound, long flag)
        !          2012: {
        !          2013:   if (flag<0 || flag>4) err(flagerr,"ideallist");
        !          2014:   return ideallistzstarall(bnf,bound,flag);
        !          2015: }
        !          2016:
        !          2017: GEN
        !          2018: ideallistzstar(GEN nf,long bound)
        !          2019: {
        !          2020:   return ideallistzstarall(nf,bound,0);
        !          2021: }
        !          2022:
        !          2023: GEN
        !          2024: ideallistzstargen(GEN nf,long bound)
        !          2025: {
        !          2026:   return ideallistzstarall(nf,bound,1);
        !          2027: }
        !          2028:
        !          2029: GEN
        !          2030: ideallistunit(GEN nf,long bound)
        !          2031: {
        !          2032:   return ideallistzstarall(nf,bound,2);
        !          2033: }
        !          2034:
        !          2035: GEN
        !          2036: ideallistunitgen(GEN nf,long bound)
        !          2037: {
        !          2038:   return ideallistzstarall(nf,bound,3);
        !          2039: }
        !          2040:
        !          2041: GEN
        !          2042: ideallist(GEN bnf,long bound)
        !          2043: {
        !          2044:   return ideallistzstarall(bnf,bound,4);
        !          2045: }
        !          2046:
        !          2047: static GEN
        !          2048: ideallist_arch(GEN nf,GEN list,GEN arch,long flun)
        !          2049: {
        !          2050:   long nba,i,j,lx,ly;
        !          2051:   GEN p1,z,p2;
        !          2052:
        !          2053:   nba=0; for (i=1; i<lg(arch); i++) if (signe(arch[i])) nba++;
        !          2054:   lx=lg(list); z=cgetg(lx,t_VEC);
        !          2055:   for (i=1; i<lx; i++)
        !          2056:   {
        !          2057:     p2=(GEN)list[i]; checkbid(p2); ly=lg(p2);
        !          2058:     p1=cgetg(ly,t_VEC); z[i]=(long)p1;
        !          2059:     for (j=1; j<ly; j++)
        !          2060:       p1[j]=(long)zidealstarinitjoinarch(nf,(GEN)p2[j],arch,nba,flun);
        !          2061:   }
        !          2062:   return z;
        !          2063: }
        !          2064:
        !          2065: static GEN
        !          2066: ideallistarchall(GEN bnf,GEN list,GEN arch,long flag)
        !          2067: {
        !          2068:   ulong av;
        !          2069:   long i,j,lp1, do_units = flag & 2;
        !          2070:   GEN nf = checknf(bnf), p1,p2,p3,racunit,funits,lu2,lu,embunit,z,y;
        !          2071:
        !          2072:   if (typ(list) != t_VEC || (do_units && lg(list) != 3))
        !          2073:     err(typeer, "ideallistarch");
        !          2074:   if (lg(list) == 1) return cgetg(1,t_VEC);
        !          2075:   if (do_units)
        !          2076:   {
        !          2077:     y = cgetg(3,t_VEC);
        !          2078:     z = (GEN)list[1];
        !          2079:     lu= (GEN)list[2]; if (typ(lu) != t_VEC) err(typeer, "ideallistarch");
        !          2080:   }
        !          2081:   else
        !          2082:   {
        !          2083:     z = list;
        !          2084:     y = lu = NULL; /* gcc -Wall */
        !          2085:   }
        !          2086:   if (typ(z) != t_VEC) err(typeer, "ideallistarch");
        !          2087:   z = ideallist_arch(nf,z,arch, flag & 1);
        !          2088:   if (!do_units) return z;
        !          2089:
        !          2090:   y[1]=(long)z; av=avma;
        !          2091:   init_units(bnf,&funits,&racunit);
        !          2092:   lu2=cgetg(lg(z),t_VEC);
        !          2093:   for (i=1; i<lg(z); i++)
        !          2094:   {
        !          2095:     p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
        !          2096:     p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
        !          2097:     for (j=1; j<lp1; j++)
        !          2098:     {
        !          2099:       embunit = logunitmatrixarch(nf,funits,racunit,(GEN)p1[j]);
        !          2100:       p3[j] = lmul(gmael(p1,j,5), vconcat((GEN)p2[j],embunit));
        !          2101:     }
        !          2102:   }
        !          2103:   y[2]=lpilecopy(av,lu2); return y;
        !          2104: }
        !          2105:
        !          2106: GEN
        !          2107: ideallistarch(GEN nf, GEN list, GEN arch)
        !          2108: {
        !          2109:   return ideallistarchall(nf,list,arch,0);
        !          2110: }
        !          2111:
        !          2112: GEN
        !          2113: ideallistarchgen(GEN nf, GEN list, GEN arch)
        !          2114: {
        !          2115:   return ideallistarchall(nf,list,arch,1);
        !          2116: }
        !          2117:
        !          2118: GEN
        !          2119: ideallistunitarch(GEN bnf,GEN list,GEN arch)
        !          2120: {
        !          2121:   return ideallistarchall(bnf,list,arch,2);
        !          2122: }
        !          2123:
        !          2124: GEN
        !          2125: ideallistunitarchgen(GEN bnf,GEN list,GEN arch)
        !          2126: {
        !          2127:   return ideallistarchall(bnf,list,arch,3);
        !          2128: }
        !          2129:
        !          2130: GEN
        !          2131: ideallistarch0(GEN nf, GEN list, GEN arch,long flag)
        !          2132: {
        !          2133:   if (!arch) arch=cgetg(1,t_VEC);
        !          2134:   if (flag<0 || flag>3) err(flagerr,"ideallistarch");
        !          2135:   return ideallistarchall(nf,list,arch,flag);
        !          2136: }

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