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

Annotation of OpenXM_contrib/pari-2.2/src/modules/subfield.c, Revision 1.1

1.1     ! noro        1: /* $Id: subfield.c,v 1.25 2001/10/01 17:19:07 bill 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: /*               SUBFIELDS OF A NUMBER FIELD                       */
        !            19: /*   J. Klueners and M. Pohst, J. Symb. Comp. (1996), vol. 11      */
        !            20: /*                                                                 */
        !            21: /*******************************************************************/
        !            22: #include "pari.h"
        !            23: extern GEN matrixpow(long n, long m, GEN y, GEN P,GEN l);
        !            24: extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
        !            25: extern GEN FpX_rand(long d1, long v, GEN p);
        !            26: extern GEN hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pe, long e);
        !            27:
        !            28: static GEN print_block_system(long N,GEN Y,long d,GEN vbs,long maxl);
        !            29:
        !            30: /* Computation of potential block systems of given size d associated to a
        !            31:  * rational prime p: give a row vector of row vectors containing the
        !            32:  * potential block systems of imprimitivity; a potential block system is a
        !            33:  * vector of row vectors (enumeration of the roots).
        !            34:  */
        !            35: #define BIL 32 /* for 64bit machines also */
        !            36: static GEN
        !            37: calc_block(long N,GEN Z,long d,GEN Y,GEN vbs,ulong maxl)
        !            38: {
        !            39:   long r,lK,i,j,t,tp,T,lpn,u,nn,lnon,lY;
        !            40:   GEN K,n,non,pn,pnon,e,Yp,Zp,Zpp;
        !            41:
        !            42:   if (DEBUGLEVEL>3)
        !            43:   {
        !            44:     long l = vbs? lg(vbs): 0;
        !            45:     fprintferr("avma = %ld, lg(Z) = %ld, lg(Y) = %ld, lg(vbs) = %ld\n",
        !            46:                avma,lg(Z),lg(Y),l);
        !            47:     if (DEBUGLEVEL > 5)
        !            48:     {
        !            49:       fprintferr("Z = %Z\n",Z);
        !            50:       fprintferr("Y = %Z\n",Y);
        !            51:       if (DEBUGLEVEL > 7) fprintferr("vbs = %Z\n",vbs);
        !            52:     }
        !            53:   }
        !            54:   r=lg(Z); lnon = min(BIL, r);
        !            55:   e    = new_chunk(BIL);
        !            56:   n    = new_chunk(r);
        !            57:   non  = new_chunk(lnon);
        !            58:   pnon = new_chunk(lnon);
        !            59:   pn   = new_chunk(lnon);
        !            60:   Zp   = cgetg(lnon,t_VEC);
        !            61:   Zpp  = cgetg(lnon,t_VEC);
        !            62:   for (i=1; i<r; i++) n[i] = lg(Z[i])-1;
        !            63:
        !            64:   K=divisors(stoi(n[1])); lK=lg(K);
        !            65:   for (i=1; i<lK; i++)
        !            66:   {
        !            67:     long  ngcd = n[1], k = itos((GEN)K[i]), dk = d*k;
        !            68:     lpn=0;
        !            69:     for (j=2; j<r; j++)
        !            70:       if (n[j]%k == 0)
        !            71:       {
        !            72:         if (++lpn >= BIL) err(talker,"overflow in calc_block");
        !            73:         pn[lpn] = n[j]; pnon[lpn] = j;
        !            74:         ngcd = cgcd(ngcd, n[j]);
        !            75:       }
        !            76:     if (dk % ngcd) continue;
        !            77:     T = 1<<lpn; if (!lpn) lpn = 1;
        !            78:     for (t=0; t<T; t++)
        !            79:     {
        !            80:       for (nn=n[1],tp=t, u=1; u<=lpn; u++,tp>>=1)
        !            81:       {
        !            82:         if (tp&1) { nn += pn[u]; e[u]=1; } else e[u]=0;
        !            83:       }
        !            84:       if (dk == nn)
        !            85:       {
        !            86:        long av=avma;
        !            87:         int Z_equal_Zp = 1;
        !            88:
        !            89:         for (j=1; j<lnon; j++) non[j]=0;
        !            90:         Zp[1]=Z[1];
        !            91:        for (u=2,j=1; j<=lpn; j++)
        !            92:          if (e[j])
        !            93:           {
        !            94:             Zp[u]=Z[pnon[j]]; non[pnon[j]]=1;
        !            95:             if (Zp[u] != Z[u]) Z_equal_Zp = 0;
        !            96:             u++;
        !            97:           }
        !            98:         setlg(Zp, u);
        !            99:         lY=lg(Y); Yp = cgetg(lY+1,t_VEC);
        !           100:         for (j=1; j<lY; j++) Yp[j]=Y[j];
        !           101:        Yp[lY]=(long)Zp;
        !           102:         if (r == u && Z_equal_Zp)
        !           103:          vbs = print_block_system(N,Yp,d,vbs,maxl);
        !           104:        else
        !           105:        {
        !           106:          for (u=1,j=2; j<r; j++)
        !           107:            if (!non[j]) Zpp[u++] = Z[j];
        !           108:           setlg(Zpp, u);
        !           109:          vbs = calc_block(N,Zpp,d,Yp,vbs,maxl);
        !           110:        }
        !           111:         if (vbs && lg(vbs) > maxl) return vbs;
        !           112:         avma=av;
        !           113:       }
        !           114:     }
        !           115:   }
        !           116:   return vbs;
        !           117: }
        !           118:
        !           119: static GEN
        !           120: potential_block_systems(long N, long d, GEN n, ulong maxl)
        !           121: {
        !           122:   long av=avma,r,i,j,k;
        !           123:   GEN p1,vbs,Z;
        !           124:
        !           125:   r=lg(n); Z=cgetg(r,t_VEC);
        !           126:   for (k=0,i=1; i<r; i++)
        !           127:   {
        !           128:     p1=cgetg(n[i]+1,t_VECSMALL); Z[i]=(long)p1;
        !           129:     for (j=1; j<=n[i]; j++) p1[j]= ++k;
        !           130:   }
        !           131:   vbs=calc_block(N,Z,d, cgetg(1,t_VEC), NULL, maxl);
        !           132:   avma=av; return vbs;
        !           133: }
        !           134:
        !           135: /* product of permutations. Put the result in perm1. */
        !           136: static void
        !           137: perm_mul(GEN perm1,GEN perm2)
        !           138: {
        !           139:   long av = avma,i, N = lg(perm1);
        !           140:   GEN perm=new_chunk(N);
        !           141:   for (i=1; i<N; i++) perm[i]=perm1[perm2[i]];
        !           142:   for (i=1; i<N; i++) perm1[i]=perm[i];
        !           143:   avma=av;
        !           144: }
        !           145:
        !           146: /* cy is a cycle; compute cy^l as a permutation */
        !           147: static GEN
        !           148: cycle_power_to_perm(GEN perm,GEN cy,long l)
        !           149: {
        !           150:   long lp,i,j,b, N = lg(perm), lcy = lg(cy)-1;
        !           151:
        !           152:   lp = l % lcy;
        !           153:   for (i=1; i<N; i++) perm[i] = i;
        !           154:   if (lp)
        !           155:   {
        !           156:     long av = avma;
        !           157:     GEN p1 = new_chunk(N);
        !           158:     b = cy[1];
        !           159:     for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
        !           160:     perm[b] = cy[1];
        !           161:     for (i=1; i<N; i++) p1[i] = perm[i];
        !           162:
        !           163:     for (j=2; j<=lp; j++) perm_mul(perm,p1);
        !           164:     avma = av;
        !           165:   }
        !           166:   return perm;
        !           167: }
        !           168:
        !           169: /* image du block system D par la permutation perm */
        !           170: static GEN
        !           171: im_block_by_perm(GEN D,GEN perm)
        !           172: {
        !           173:   long i,j,lb,lcy;
        !           174:   GEN Dn,cy,p1;
        !           175:
        !           176:   lb=lg(D); Dn=cgetg(lb,t_VEC);
        !           177:   for (i=1; i<lb; i++)
        !           178:   {
        !           179:     cy=(GEN)D[i]; lcy=lg(cy);
        !           180:     Dn[i]=lgetg(lcy,t_VECSMALL); p1=(GEN)Dn[i];
        !           181:     for (j=1; j<lcy; j++) p1[j] = perm[cy[j]];
        !           182:   }
        !           183:   return Dn;
        !           184: }
        !           185:
        !           186: /* cy is a cycle; recturn cy(a) */
        !           187: static long
        !           188: im_by_cy(long a,GEN cy)
        !           189: {
        !           190:   long k, l = lg(cy);
        !           191:
        !           192:   k=1; while (k<l && cy[k] != a) k++;
        !           193:   if (k == l) return a;
        !           194:   k++; if (k == l) k = 1;
        !           195:   return cy[k];
        !           196: }
        !           197:
        !           198: /* vbs[0] = current cardinality+1, vbs[1] = max number of elts */
        !           199: static GEN
        !           200: append_vbs(GEN vbs, GEN D)
        !           201: {
        !           202:   long l,maxl,i,j,n, lD = lg(D);
        !           203:   GEN Dn, last;
        !           204:
        !           205:   n = 0; for (i=1; i<lD; i++) n += lg(D[i]);
        !           206:   Dn = (GEN)gpmalloc((lD + n) * sizeof(long));
        !           207:   last = Dn + lD; Dn[0] = D[0];
        !           208:   for (i=1; i<lD; i++)
        !           209:   {
        !           210:     GEN cy = (GEN)D[i], cn = last;
        !           211:     for (j=0; j<lg(cy); j++) cn[j] = cy[j];
        !           212:     Dn[i] = (long)cn; last = cn + j;
        !           213:   }
        !           214:
        !           215:   if (!vbs)
        !           216:   {
        !           217:     maxl = 1024;
        !           218:     vbs = (GEN)gpmalloc((2 + maxl)*sizeof(GEN));
        !           219:     *vbs = maxl; vbs++; setlg(vbs, 1);
        !           220:   }
        !           221:   l = lg(vbs); maxl = vbs[-1];
        !           222:   if (l == maxl)
        !           223:   {
        !           224:     vbs = (GEN)gprealloc((void*)(vbs-1), (2 + (maxl<<1))*sizeof(GEN),
        !           225:                                          (2 + maxl)*sizeof(GEN));
        !           226:     *vbs = maxl<<1; vbs++; setlg(vbs, 1);
        !           227:   }
        !           228:   if (DEBUGLEVEL>5) fprintferr("appending D = %Z\n",D);
        !           229:   vbs[l] = (long)Dn; setlg(vbs, l+1); return vbs;
        !           230: }
        !           231:
        !           232: GEN
        !           233: myconcat(GEN D, long a)
        !           234: {
        !           235:   long i,l = lg(D);
        !           236:   GEN x = cgetg(l+1,t_VECSMALL);
        !           237:   for (i=1; i<l; i++) x[i]=D[i];
        !           238:   x[l] = a; return x;
        !           239: }
        !           240:
        !           241: void
        !           242: myconcat2(GEN D, GEN a)
        !           243: {
        !           244:   long i,l = lg(D), m = lg(a);
        !           245:   GEN x = D + (l-1);
        !           246:   for (i=1; i<m; i++) x[i]=a[i];
        !           247:   setlg(D, l+m-1);
        !           248: }
        !           249:
        !           250: static GEN
        !           251: print_block_system(long N,GEN Y,long d,GEN vbs,long maxl)
        !           252: {
        !           253:   long a,i,j,l,ll,*k,*n,lp,**e,u,v,*t,ns, r = lg(Y);
        !           254:   GEN D,De,Z,cyperm,perm,p1,empty;
        !           255:
        !           256:   if (DEBUGLEVEL>5) fprintferr("Y = %Z\n",Y);
        !           257:   empty = cgetg(1,t_VEC);
        !           258:   n = new_chunk(N+1);
        !           259:   D = cgetg(N+1, t_VEC); setlg(D,1);
        !           260:   t=new_chunk(r+1); k=new_chunk(r+1); Z=cgetg(r+1,t_VEC);
        !           261:   for (ns=0,i=1; i<r; i++)
        !           262:   {
        !           263:     GEN Yi = (GEN)Y[i], cy;
        !           264:     long ki = 0, si = lg(Yi)-1;
        !           265:
        !           266:     for (j=1; j<=si; j++) { n[j]=lg(Yi[j])-1; ki += n[j]; }
        !           267:     ki /= d;
        !           268:     De=cgetg(ki+1,t_VEC);
        !           269:     for (j=1; j<=ki; j++) De[j]=(long)empty;
        !           270:     for (j=1; j<=si; j++)
        !           271:     {
        !           272:       cy = (GEN)Yi[j]; a = cy[1];
        !           273:       for (l=1,lp=0; l<=n[j]; l++)
        !           274:       {
        !           275:         lp++; if (lp>ki) lp = 1;
        !           276:         a = im_by_cy(a, cy);
        !           277:         De[lp] = (long)myconcat((GEN)De[lp], a);
        !           278:       }
        !           279:     }
        !           280:     if (si>1 && ki>1)
        !           281:     {
        !           282:       ns++; t[ns]=si-1; k[ns]=ki;
        !           283:       Z[ns]=lgetg(si,t_VEC); p1=(GEN)Z[ns];
        !           284:       for (j=2; j<=si; j++) p1[j-1]=Yi[j];
        !           285:     }
        !           286:     myconcat2(D,De);
        !           287:   }
        !           288:   if (DEBUGLEVEL>2) { fprintferr("\nns = %ld\n",ns); flusherr(); }
        !           289:   if (!ns) return append_vbs(vbs,D);
        !           290:
        !           291:   setlg(Z, ns+1);
        !           292:   e=(long**)new_chunk(ns+1);
        !           293:   for (i=1; i<=ns; i++)
        !           294:   {
        !           295:     e[i]=new_chunk(t[i]+1);
        !           296:     for (j=1; j<=t[i]; j++) e[i][j]=0;
        !           297:   }
        !           298:   cyperm = cgetg(N+1,t_VEC);
        !           299:   perm = cgetg(N+1,t_VEC); i=ns;
        !           300:   do
        !           301:   {
        !           302:     long av = avma;
        !           303:     if (DEBUGLEVEL>5)
        !           304:     {
        !           305:       for (l=1; l<=ns; l++)
        !           306:       {
        !           307:        for (ll=1; ll<=t[l]; ll++)
        !           308:          fprintferr("e[%ld][%ld] = %ld, ",l,ll,e[l][ll]);
        !           309:        fprintferr("\n");
        !           310:       }
        !           311:       fprintferr("\n"); flusherr();
        !           312:     }
        !           313:     for (u=1; u<=N; u++) perm[u]=u;
        !           314:     for (u=1; u<=ns; u++)
        !           315:       for (v=1; v<=t[u]; v++)
        !           316:        perm_mul(perm, cycle_power_to_perm(cyperm,gmael(Z,u,v),e[u][v]));
        !           317:     vbs = append_vbs(vbs, im_block_by_perm(D,perm));
        !           318:     if (lg(vbs) > maxl) return vbs;
        !           319:     avma = av;
        !           320:
        !           321:     e[ns][t[ns]]++;
        !           322:     if (e[ns][t[ns]] >= k[ns])
        !           323:     {
        !           324:       j=t[ns]-1;
        !           325:       while (j>=1 && e[ns][j] == k[ns]-1) j--;
        !           326:       if (j>=1) { e[ns][j]++; for (l=j+1; l<=t[ns]; l++) e[ns][l]=0; }
        !           327:       else
        !           328:       {
        !           329:        i=ns-1;
        !           330:        while (i>=1)
        !           331:        {
        !           332:          j=t[i];
        !           333:          while (j>=1 && e[i][j] == k[i]-1) j--;
        !           334:          if (j<1) i--;
        !           335:           else
        !           336:          {
        !           337:            e[i][j]++;
        !           338:            for (l=j+1; l<=t[i]; l++) e[i][l]=0;
        !           339:            for (ll=i+1; ll<=ns; ll++)
        !           340:               for (l=1; l<=t[ll]; l++) e[ll][l]=0;
        !           341:             break;
        !           342:          }
        !           343:        }
        !           344:       }
        !           345:     }
        !           346:   }
        !           347:   while (i>0);
        !           348:   return vbs;
        !           349: }
        !           350:
        !           351: static GEN
        !           352: polsimplify(GEN x)
        !           353: {
        !           354:   long i,lx = lgef(x);
        !           355:   for (i=2; i<lx; i++)
        !           356:     if (typ(x[i]) == t_POL) x[i] = (long)constant_term(x[i]);
        !           357:   return x;
        !           358: }
        !           359:
        !           360: /* return 0 if |g[i]| > M[i] for some i; 1 otherwise */
        !           361: static long
        !           362: ok_coeffs(GEN g,GEN M)
        !           363: {
        !           364:   long i, lg = lgef(g)-1; /* g is monic, and cst term is ok */
        !           365:   for (i=3; i<lg; i++)
        !           366:     if (absi_cmp((GEN)g[i], (GEN)M[i]) > 0) return 0;
        !           367:   return 1;
        !           368: }
        !           369:
        !           370: /* Return a polynomial g defining a potential subfield, or
        !           371:  * 0: if p | d(g)
        !           372:  * 1: if coeffs of g are too large
        !           373:  * 2: same, detected by cheap d-1 test */
        !           374: static GEN
        !           375: cand_for_subfields(GEN A,GEN DATA,GEN *ptlistdelta)
        !           376: {
        !           377:   long N,m,i,j,d,lf;
        !           378:   GEN M,T,pe,p,pol,fhk,g;
        !           379:   GEN _d_1_term,delta,listdelta,whichdelta,firstroot;
        !           380:
        !           381:   pol=(GEN)DATA[1];
        !           382:   p = (GEN)DATA[2];
        !           383:   pe= (GEN)DATA[3];
        !           384:   T = (GEN)DATA[4];
        !           385:   fhk=(GEN)DATA[5];
        !           386:   M = (GEN)DATA[8]; N=degpol(pol); m=lg(A)-1; d=N/m; /* m | M */
        !           387:   firstroot = (GEN)DATA[13];
        !           388:
        !           389:   delta = cgetg(m+1,t_VEC);
        !           390:   lf = lg(firstroot);
        !           391:   listdelta = cgetg(lf, t_VEC);
        !           392:   whichdelta = cgetg(N+1, t_VECSMALL);
        !           393:   _d_1_term = gzero;
        !           394:   for (i=1; i<=m; i++)
        !           395:   {
        !           396:     GEN Ai = (GEN)A[i], p1 = gun;
        !           397:     for (j=1; j<=d; j++)
        !           398:       p1 = FpXQ_mul(p1, (GEN)fhk[Ai[j]], T,pe);
        !           399:     delta[i] = (long)p1;
        !           400:     if (DEBUGLEVEL>2) fprintferr("delta[%ld] = %Z\n",i,p1);
        !           401:     /* g = prod (X - delta[i])
        !           402:      * if g o h = 0 (pol), we'll have h(Ai[j]) = delta[i] for all j */
        !           403:     /* fk[k] belongs to block number whichdelta[k] */
        !           404:     for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
        !           405:     if (typ(p1) == t_POL) p1 = constant_term(p1);
        !           406:     _d_1_term = addii(_d_1_term, p1);
        !           407:   }
        !           408:   _d_1_term = centermod(_d_1_term, pe); /* Tr(g) */
        !           409:   if (absi_cmp(_d_1_term, (GEN)M[3]) > 0) return gdeux; /* d-1 test */
        !           410:   g = FqV_roots_to_pol(delta, T, pe, 0);
        !           411:   g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
        !           412:   if (DEBUGLEVEL>2) fprintferr("pol. found = %Z\n",g);
        !           413:   if (!ok_coeffs(g,M)) return gun;
        !           414:   if (!FpX_is_squarefree(g, p)) return gzero;
        !           415:
        !           416:   for (i=1; i<lf; i++) listdelta[i] = delta[whichdelta[firstroot[i]]];
        !           417:   *ptlistdelta = listdelta; return g;
        !           418: }
        !           419:
        !           420: /* return U list of polynomials s.t U[i] = 1 mod fk[i] and 0 mod fk[j] for all
        !           421:  * other j */
        !           422: static GEN
        !           423: get_bezout(GEN pol, GEN fk, GEN p)
        !           424: {
        !           425:   GEN A,B,d,u,v,bezout_coeff;
        !           426:   long i, l = lg(fk);
        !           427:
        !           428:   pol = FpX_red(pol, p);
        !           429:   bezout_coeff = cgetg(l, t_VEC);
        !           430:   for (i=1; i<l; i++)
        !           431:   {
        !           432:     A = (GEN)fk[i];
        !           433:     B = FpX_div(pol, A, p);
        !           434:     d = FpX_extgcd(A,B,p, &u, &v);
        !           435:     if (degpol(d) > 0) err(talker, "relatively prime polynomials expected");
        !           436:     d = constant_term(d);
        !           437:     if (!gcmp1(d))
        !           438:     {
        !           439:       d = mpinvmod(d, p);
        !           440:       v = FpX_Fp_mul(v,d, p);
        !           441:     }
        !           442:     bezout_coeff[i] = (long)FpX_mul(B,v, p);
        !           443:   }
        !           444:   return bezout_coeff;
        !           445: }
        !           446:
        !           447: /* assume x in Fq, return Tr_{Fq/Fp}(x) as a t_INT */
        !           448: static GEN
        !           449: trace(GEN x, GEN Trq, GEN p)
        !           450: {
        !           451:   long i, l;
        !           452:   GEN s;
        !           453:   if (typ(x) == t_INT) return modii(mulii(x, (GEN)Trq[1]), p);
        !           454:   l = lgef(x)-1; if (l == 1) return gzero;
        !           455:   x++; s = mulii((GEN)x[1], (GEN)Trq[1]);
        !           456:   for (i=2; i<l; i++)
        !           457:     s = addii(s, mulii((GEN)x[i], (GEN)Trq[i]));
        !           458:   return modii(s, p);
        !           459: }
        !           460:
        !           461: /* assume x in Fq[X], return Tr_{Fq[X]/Fp[X]}(x), varn(X) = 0 */
        !           462: static GEN
        !           463: poltrace(GEN x, GEN Trq, GEN p)
        !           464: {
        !           465:   long i,l;
        !           466:   GEN y;
        !           467:   if (typ(x) == t_INT || varn(x) != 0) return trace(x, Trq, p);
        !           468:   l = lgef(x); y = cgetg(l,t_POL); y[1]=x[1];
        !           469:   for (i=2; i<l; i++) y[i] = (long)trace((GEN)x[i],Trq,p);
        !           470:   return y;
        !           471: }
        !           472:
        !           473: /* Find h in Fp[X] such that h(a[i]) = listdelta[i] for all modular factors
        !           474:  * ff[i], where a[i] is a fixed root of ff[i] in Fq = Z[Y]/(p,T) [namely the
        !           475:  * first one in Fp_factor_irred output]. Let f = ff[i], A the given root, then
        !           476:  * h mod f is Tr_Fq/Fp ( h(A) f(X)/(X-A)f'(A) ), most of the expression being
        !           477:  * precomputed. The complete h is recovered via chinese remaindering */
        !           478: static GEN
        !           479: chinese_retrieve_pol(GEN DATA, GEN listdelta)
        !           480: {
        !           481:   GEN interp,Trk,bezoutC,T,p, S,p1;
        !           482:   long i,l;
        !           483:   p = (GEN)DATA[2];
        !           484:   T = (GEN)DATA[4];
        !           485:   interp = (GEN)DATA[10];
        !           486:   Trk = (GEN)DATA[11];
        !           487:   bezoutC = (GEN)DATA[12];
        !           488:
        !           489:   S = NULL; l = lg(interp);
        !           490:   for (i=1; i<l; i++)
        !           491:   { /* h(firstroot[i]) = listdelta[i] */
        !           492:     p1 = FpXQX_FpXQ_mul((GEN)interp[i], (GEN)listdelta[i], T,p);
        !           493:     p1 = poltrace(p1, (GEN)Trk[i], p);
        !           494:     p1 = gmul(p1, (GEN)bezoutC[i]);
        !           495:     S = S? gadd(S,p1): p1;
        !           496:   }
        !           497:   return FpX_res(FpX_red(S, p), FpX_red((GEN)DATA[1],p), p);
        !           498: }
        !           499:
        !           500: /* g in Z[X] potentially defines a subfield of Q[X]/f. It is a subfield iff A
        !           501:  * (cf cand_for_subfields) was a block system; then there
        !           502:  * exists h in Q[X] such that f | g o h. listdelta determines h s.t f | g o h
        !           503:  * in Fp[X] (cf chinese_retrieve_pol). We try to lift it. */
        !           504: static GEN
        !           505: embedding_of_potential_subfields(GEN g,GEN DATA,GEN listdelta)
        !           506: {
        !           507:   GEN TR,w0_Q,w0,w1_Q,w1,wpow,h0,gp,T,q2,q,p,ind,maxp,a;
        !           508:   long rt;
        !           509:   ulong av;
        !           510:
        !           511:   T = (GEN)DATA[1]; rt = brent_kung_optpow(degpol(T), 2);
        !           512:   p = (GEN)DATA[2];
        !           513:   maxp = (GEN)DATA[7];
        !           514:   ind = (GEN)DATA[9];
        !           515:   gp = derivpol(g); av = avma;
        !           516:   w0 = chinese_retrieve_pol(DATA,listdelta);
        !           517:   w0_Q = centermod(gmul(w0,ind), p);
        !           518:   h0 = FpXQ_inv(FpX_FpXQ_compo(gp,w0, T,p), T,p); /* = 1/g'(w0) mod (T,p) */
        !           519:   wpow = NULL; q = sqri(p);
        !           520:   for(;;)
        !           521:   {/* Given g,w0,h0 in Z[x], s.t. h0.g'(w0) = 1 and g(w0) = 0 mod (T,p), find
        !           522:     * [w1,h1] satisfying the same conditions mod p^2, [w1,h1] = [w0,h0] (mod p)
        !           523:     * (cf. Dixon: J. Austral. Math. Soc., Series A, vol.49, 1990, p.445) */
        !           524:     if (DEBUGLEVEL>1) fprintferr("lifting embedding mod p = %Z\n",q);
        !           525:
        !           526:     /* w1 := w0 - h0 g(w0) mod (T,q) */
        !           527:     if (wpow)
        !           528:       a = FpX_FpXQV_compo(g,wpow, T,q);
        !           529:     else
        !           530:       a = FpX_FpXQ_compo(g,w0, T,q); /* first time */
        !           531:     /* now, a = 0 (p) */
        !           532:     a = gmul(gneg(h0), gdivexact(a, p));
        !           533:     w1 = gadd(w0, gmul(p, FpX_res(a, T,p)));
        !           534:
        !           535:     w1_Q = centermod(gmul(w1, resii(ind,q)), q);
        !           536:     if (gegal(w1_Q, w0_Q) || cmpii(q,maxp) > 0)
        !           537:     {
        !           538:       GEN G = gcmp1(ind)? g: ZX_rescale_pol(g,ind);
        !           539:       if (gcmp0(poleval(G, gmodulcp(w1_Q,T)))) break;
        !           540:     }
        !           541:     if (cmpii(q, maxp) > 0)
        !           542:     {
        !           543:       if (DEBUGLEVEL) fprintferr("coeff too big for embedding\n");
        !           544:       return NULL;
        !           545:     }
        !           546:     {
        !           547:       GEN *gptr[5]; gptr[0]=&w1; gptr[1]=&h0; gptr[2]=&w1_Q;
        !           548:       gptr[3]=&q; gptr[4]=&p;
        !           549:       gerepilemany(av,gptr,5);
        !           550:     }
        !           551:
        !           552:     q2 = sqri(q);
        !           553:     wpow = FpXQ_powers(w1, rt, T, q2);
        !           554:     /* h0 := h0 * (2 - h0 g'(w1)) mod (T,q)
        !           555:      *     = h0 + h0 * (1 - h0 g'(w1)) */
        !           556:     a = gmul(gneg(h0), FpX_FpXQV_compo(gp, FpXV_red(wpow,q),T,q));
        !           557:     a = ZX_s_add(FpX_res(a, T,q), 1); /* 1 - h0 g'(w1) = 0 (p) */
        !           558:     a = gmul(h0, gdivexact(a, p));
        !           559:     h0 = gadd(h0, gmul(p, FpX_res(a, T,p)));
        !           560:     w0 = w1; w0_Q = w1_Q; p = q; q = q2;
        !           561:   }
        !           562:   TR = (GEN)DATA[14];
        !           563:   if (!gcmp0(TR)) w1_Q = poleval(w1_Q, gadd(polx[0], TR));
        !           564:   return gdiv(w1_Q,ind);
        !           565: }
        !           566:
        !           567: static GEN
        !           568: choose_prime(GEN pol,GEN dpol,long d,GEN *ptff,GEN *ptlistpotbl, long *ptlcm)
        !           569: {
        !           570:   ulong av;
        !           571:   long j,k,oldllist,llist,r,lcm,oldlcm,N,pp;
        !           572:   GEN p,listpotbl,oldlistpotbl,ff,oldff,n,oldn;
        !           573:   byteptr di=diffptr;
        !           574:
        !           575:   if (DEBUGLEVEL) timer2();
        !           576:   di++; p = stoi(2); N = degpol(pol);
        !           577:   while (p[2]<=N) p[2] += *di++;
        !           578:   oldllist = 100000;
        !           579:   oldlcm = 0;
        !           580:   oldlistpotbl = oldff = oldn = NULL; pp = 0; /* gcc -Wall */
        !           581:   av = avma;
        !           582:   for(k=1; k<11 || !oldn; k++,p[2]+= *di++,avma = av)
        !           583:   {
        !           584:     while (!smodis(dpol,p[2])) p[2] += *di++;
        !           585:     if (k > 50) err(talker,"sorry, too many block systems in nfsubfields");
        !           586:     ff=(GEN)factmod(pol,p)[1]; r=lg(ff)-1;
        !           587:     if (r == 1 || r == N) continue;
        !           588:
        !           589:     n = cgetg(r+1, t_VECSMALL);
        !           590:     lcm = n[1] = degpol(ff[1]);
        !           591:     for (j=2; j<=r; j++) { n[j] = degpol(ff[j]); lcm = clcm(lcm,n[j]); }
        !           592:     if (lcm < oldlcm) continue; /* false when oldlcm = 0 */
        !           593:     if (r >= BIL) { err(warner,"subfields: overflow in calc_block"); continue; }
        !           594:     if (DEBUGLEVEL) fprintferr("p = %ld,\tlcm = %ld,\torbits: %Z\n",p[2],lcm,n);
        !           595:     if (oldn && gegal(n,oldn)) continue;
        !           596:
        !           597:     listpotbl = potential_block_systems(N,d,n, oldllist);
        !           598:     if (!listpotbl) { oldlistpotbl = NULL; pp = p[2]; break; }
        !           599:     llist = lg(listpotbl)-1;
        !           600:     if (llist >= oldllist)
        !           601:     {
        !           602:       if (DEBUGLEVEL) msgtimer("#pbs >= %ld [aborted]",oldllist);
        !           603:       for (j=1; j<llist; j++) free((void*)listpotbl[j]);
        !           604:       free((void*)(listpotbl-1)); continue;
        !           605:     }
        !           606:     oldllist = llist; oldlistpotbl = listpotbl;
        !           607:     pp = p[2]; oldff = ff; oldlcm = lcm; oldn = n;
        !           608:     if (DEBUGLEVEL) msgtimer("#pbs = %ld",llist);
        !           609:     av = avma;
        !           610:   }
        !           611:   if (DEBUGLEVEL)
        !           612:   {
        !           613:     fprintferr("Chosen prime: p = %ld\n",pp);
        !           614:     if (DEBUGLEVEL>2)
        !           615:       fprintferr("Potential block systems of size %ld: %Z\n", d,oldlistpotbl);
        !           616:     flusherr();
        !           617:   }
        !           618:   if (oldff) oldff = lift_intern(oldff);
        !           619:   *ptlistpotbl=oldlistpotbl; *ptff=oldff; *ptlcm=oldlcm; return stoi(pp);
        !           620: }
        !           621:
        !           622: static GEN
        !           623: bound_for_coeff(long m,GEN rr, GEN *maxroot)
        !           624: {
        !           625:   long i,r1, lrr=lg(rr);
        !           626:   GEN p1,b1,b2,B,M, C = matpascal(m-1);
        !           627:
        !           628:   for (r1=0; typ(rr[r1+1]) == t_REAL; r1++) /* empty */;
        !           629:
        !           630:   rr = gabs(rr,0); *maxroot = vecmax(rr);
        !           631:   for (i=1; i<lrr; i++)
        !           632:     if (gcmp((GEN)rr[i], gun) < 0) rr[i] = un;
        !           633:   for (b1=gun,i=1; i<=r1; i++) b1 = gmul(b1, (GEN)rr[i]);
        !           634:   for (b2=gun    ; i<lrr; i++) b2 = gmul(b2, (GEN)rr[i]);
        !           635:   B = gmul(b1, gsqr(b2)); /* Mahler measure */
        !           636:   M = cgetg(m+2, t_VEC); M[1]=M[2]=zero; /* unused */
        !           637:   for (i=1; i<m; i++)
        !           638:   {
        !           639:     p1 = gadd(gmul(gcoeff(C, m, i+1), B),/* binom(m-1, i)   */
        !           640:               gcoeff(C, m, i));          /* binom(m-1, i-1) */
        !           641:     M[i+2] = (long)ceil_safe(p1);
        !           642:   }
        !           643:   return M;
        !           644: }
        !           645:
        !           646: static GEN
        !           647: init_traces(GEN ff, GEN T, GEN p)
        !           648: {
        !           649:   long N = degpol(T),i,j,k, r = lg(ff);
        !           650:   GEN Frob = matrixpow(N,N, FpXQ_pow(polx[varn(T)],p, T,p), T,p);
        !           651:   GEN y,p1,p2,Trk,pow,pow1;
        !           652:
        !           653:   k = degpol(ff[r-1]); /* largest degree in modular factorization */
        !           654:   pow = cgetg(k+1, t_VEC);
        !           655:   pow[1] = (long)zero; /* dummy */
        !           656:   pow[2] = (long)Frob;
        !           657:   pow1= cgetg(k+1, t_VEC); /* 1st line */
        !           658:   for (i=3; i<=k; i++)
        !           659:     pow[i] = (long)FpM_mul((GEN)pow[i-1], Frob, p);
        !           660:   pow1[1] = (long)zero; /* dummy */
        !           661:   for (i=2; i<=k; i++)
        !           662:   {
        !           663:     p1 = cgetg(N+1, t_VEC);
        !           664:     pow1[i] = (long)p1; p2 = (GEN)pow[i];
        !           665:     for (j=1; j<=N; j++) p1[j] = coeff(p2,1,j);
        !           666:   }
        !           667:   p1 = cgetg(N+1,t_VEC); p1[1] = un;
        !           668:   for (i=2; i<=N; i++) p1[i] = zero;
        !           669:   /* Trk[i] = line 1 of x -> x + x^p + ... + x^{p^(i-1)} */
        !           670:   Trk = pow; /* re-use (destroy) pow */
        !           671:   Trk[1] = (long)p1;
        !           672:   for (i=2; i<=k; i++)
        !           673:     Trk[i] = ladd((GEN)Trk[i-1], (GEN)pow1[i]);
        !           674:   y = cgetg(r, t_VEC);
        !           675:   for (i=1; i<r; i++) y[i] = Trk[degpol(ff[i])];
        !           676:   return y;
        !           677: }
        !           678:
        !           679: /* return C in Z[X]/(p,T), C[ D[1] ] = 1, C[ D[i] ] = 0 otherwise. H is the
        !           680:  * list of degree 1 polynomials X - D[i]  (come directly from factorization) */
        !           681: static GEN
        !           682: interpol(GEN H, GEN T, GEN p)
        !           683: {
        !           684:   long i, m = lg(H);
        !           685:   GEN X = polx[0],d,p1,p2,a;
        !           686:
        !           687:   p1=polun[0]; p2=gun; a = gneg(constant_term(H[1])); /* = D[1] */
        !           688:   for (i=2; i<m; i++)
        !           689:   {
        !           690:     d = constant_term(H[i]); /* -D[i] */
        !           691:     p1 = FpXQX_mul(p1,gadd(X,d), T,p);
        !           692:     p2 = FpXQ_mul(p2, gadd(a,d), T,p);
        !           693:   }
        !           694:   p2 = FpXQ_inv(p2, T,p);
        !           695:   return FpXQX_FpXQ_mul(p1,p2, T,p);
        !           696: }
        !           697:
        !           698: static GEN
        !           699: roots_from_deg1(GEN x)
        !           700: {
        !           701:   long i,l = lg(x);
        !           702:   GEN r = cgetg(l,t_VEC);
        !           703:   for (i=1; i<l; i++) r[i] = lneg(constant_term(x[i]));
        !           704:   return r;
        !           705: }
        !           706:
        !           707: struct poldata
        !           708: {
        !           709:   GEN pol;
        !           710:   GEN dis; /* |disc(pol)| */
        !           711:   GEN roo; /* roots(pol) */
        !           712:   GEN den; /* multiple of index(pol) */
        !           713: };
        !           714:
        !           715: /* ff = factmod(nf[1], p), d = requested degree for subfield. Return DATA,
        !           716:  * valid for given nf, p and d
        !           717:  * 1: polynomial nf[1],
        !           718:  * 2: prime p,
        !           719:  * 3: exponent e (for Hensel lifts) such that p^e > max(M),
        !           720:  * 4: polynomial defining the field F_(p^nn),
        !           721:  * 5: Hensel lift to precision p^e  of DATA[6]
        !           722:  * 6: roots of nf[1] in F_(p^nn),
        !           723:  * 7: Hadamard bound for coefficients of h(x) such that g o h = 0 mod nf[1].
        !           724:  * 8: bound M for polynomials defining subfields x DATA[9]
        !           725:  * 9: multiple of index of nf
        !           726:  *10: *[i] = interpolation polynomial for ff[i] [= 1 on the first root
        !           727:       firstroot[i], 0 on the others]
        !           728:  *11: Trk used to compute traces (cf poltrace)
        !           729:  *12: Bezout coefficients associated to the ff[i]
        !           730:  *13: *[i] = index of first root of ff[i] (in DATA[6])
        !           731:  *14: number of polynomial changes (translations) */
        !           732: static GEN
        !           733: compute_data(GEN DATA, struct poldata PD, long d, GEN ff, GEN T,GEN p)
        !           734: {
        !           735:   long i,j,l,e,N;
        !           736:   GEN den,roo,pe,p1,p2,fk,fhk,MM,maxroot,pol,interp,bezoutC;
        !           737:
        !           738:   if (DEBUGLEVEL>1) { fprintferr("Entering compute_data()\n\n"); flusherr(); }
        !           739:   pol = PD.pol; N = degpol(pol);
        !           740:   roo = PD.roo;
        !           741:   den = PD.den;
        !           742:   if (DATA) /* update (translate) an existing DATA */
        !           743:   {
        !           744:     GEN Xm1 = gsub(polx[varn(pol)], gun);
        !           745:     GEN TR = addis((GEN)DATA[14],1);
        !           746:     DATA[14] = (long)TR;
        !           747:     pol = poleval(pol, gsub(polx[varn(pol)], TR));
        !           748:     p1 = dummycopy(roo); l = lg(p1);
        !           749:     for (i=1; i<l; i++) p1[i] = ladd(TR, (GEN)p1[i]);
        !           750:     roo = p1;
        !           751:
        !           752:     fk = (GEN)DATA[6]; l=lg(fk);
        !           753:     for (i=1; i<l; i++) fk[i] = lsub(Xm1, (GEN)fk[i]);
        !           754:
        !           755:     interp = (GEN)DATA[10];
        !           756:     bezoutC = (GEN)DATA[12]; l = lg(interp);
        !           757:     for (i=1; i<l; i++)
        !           758:     {
        !           759:       if (degpol(interp[i]) > 0) /* do not turn polun[0] into gun */
        !           760:       {
        !           761:         p1 = poleval((GEN)interp[i], Xm1);
        !           762:         interp[i] = (long)FpXQX_red(p1, NULL,p);
        !           763:       }
        !           764:       if (degpol(bezoutC[i]) > 0)
        !           765:       {
        !           766:         p1 = poleval((GEN)bezoutC[i], Xm1);
        !           767:         bezoutC[i] = (long)FpXQX_red(p1, NULL,p);
        !           768:       }
        !           769:     }
        !           770:   }
        !           771:   else
        !           772:   {
        !           773:     GEN firstroot;
        !           774:     long r = lg(ff);
        !           775:     DATA = cgetg(15,t_VEC);
        !           776:     DATA[2] = (long)p;
        !           777:     DATA[4] = (long)T;
        !           778:     interp = cgetg(r, t_VEC);
        !           779:     firstroot = cgetg(r, t_VECSMALL);
        !           780:     fk = cgetg(N+1,t_VEC);
        !           781:     for (l=1,j=1; j<r; j++)
        !           782:     { /* compute roots and fix ordering (Frobenius cycles) */
        !           783:       p1 = Fp_factor_irred((GEN)ff[j], p, (GEN)DATA[4]);
        !           784:       interp[j] = (long)interpol(p1,T,p);
        !           785:       firstroot[j] = l;
        !           786:       for (i=1; i<lg(p1); i++,l++) fk[l] = p1[i];
        !           787:     }
        !           788:     DATA[9] = (long)PD.den;
        !           789:     DATA[10]= (long)interp;
        !           790:     DATA[11]= (long)init_traces(ff, T,p);
        !           791:     DATA[12]= (long)get_bezout(pol,ff,p);
        !           792:     DATA[13]= (long)firstroot;
        !           793:     DATA[14]= zero;
        !           794:   }
        !           795:   DATA[1] = (long)pol;
        !           796:   MM = bound_for_coeff(d, roo, &maxroot);
        !           797:   MM = gmul2n(MM,1);
        !           798:   DATA[8] = (long)MM;
        !           799:   e = logint(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [for d-1 test] */
        !           800:   DATA[3] = (long)pe;
        !           801:   DATA[6] = (long)roots_from_deg1(fk);
        !           802:   fhk = hensel_lift_fact(pol,fk,(GEN)DATA[4],p,pe,e);
        !           803:   DATA[5] = (long)roots_from_deg1(fhk);
        !           804:
        !           805:   p1 = gmul(stoi(N), gsqrt(gpuigs(stoi(N-1),N-1),DEFAULTPREC));
        !           806:   p2 = gpuigs(maxroot, N/d + N*(N-1)/2);
        !           807:   p1 = gdiv(gmul(p1,p2), gsqrt(PD.dis,DEFAULTPREC));
        !           808:   p1 = shifti(ceil_safe(p1), 1);
        !           809:   DATA[7] = lmulii(p1, PD.den);
        !           810:
        !           811:   if (DEBUGLEVEL>1) {
        !           812:     fprintferr("f = %Z\n",DATA[1]);
        !           813:     fprintferr("p = %Z, lift to p^%ld\n",DATA[2], e);
        !           814:     fprintferr("Fq defined by %Z\n",DATA[4]);
        !           815:     fprintferr("2 * Hadamard bound * ind = %Z\n",DATA[7]);
        !           816:     fprintferr("2 * M = %Z\n",DATA[8]);
        !           817:   }
        !           818:   return DATA;
        !           819: }
        !           820:
        !           821: /* g = polynomial, h = embedding. Return [[g,h]] */
        !           822: static GEN
        !           823: _subfield(GEN g, GEN h)
        !           824: {
        !           825:   GEN x = cgetg(3,t_VEC);
        !           826:   x[1] = (long)g;
        !           827:   x[2] = (long)h; return _vec(x);
        !           828: }
        !           829:
        !           830: /* subfields of degree d */
        !           831: static GEN
        !           832: subfields_of_given_degree(struct poldata PD,long d)
        !           833: {
        !           834:   ulong av,av2;
        !           835:   long llist,i,nn;
        !           836:   GEN listpotbl,ff,A,CSF,ESF,LSB,p,T,DATA,listdelta;
        !           837:   GEN pol = PD.pol, dpol = PD.dis;
        !           838:
        !           839:   if (DEBUGLEVEL) fprintferr("\n*** Look for subfields of degree %ld\n\n", d);
        !           840:   av = avma;
        !           841:   p = choose_prime(pol,dpol,degpol(pol)/d,&ff,&listpotbl,&nn);
        !           842:   if (!listpotbl) { avma=av; return cgetg(1,t_VEC); }
        !           843:   T = lift_intern(ffinit(p,nn, fetch_var()));
        !           844:   DATA = NULL; LSB = cgetg(1,t_VEC);
        !           845:   i = 1; llist = lg(listpotbl);
        !           846: CHANGE:
        !           847:   DATA = compute_data(DATA,PD,d, ff,T,p);
        !           848:   for (; i<llist; i++)
        !           849:   {
        !           850:     av2 = avma; A = (GEN)listpotbl[i];
        !           851:     if (DEBUGLEVEL > 1) fprintferr("\n* Potential block # %ld: %Z\n",i,A);
        !           852:     CSF = cand_for_subfields(A,DATA,&listdelta);
        !           853:     if (typ(CSF)==t_INT)
        !           854:     {
        !           855:       avma = av2;
        !           856:       if (DEBUGLEVEL > 1) switch(itos(CSF))
        !           857:       {
        !           858:         case 0: fprintferr("changing f(x): p divides disc(g(x))\n"); break;
        !           859:         case 1: fprintferr("coeff too big for pol g(x)\n"); break;
        !           860:         case 2: fprintferr("d-1 test failed\n"); break;
        !           861:       }
        !           862:       if (CSF==gzero) goto CHANGE;
        !           863:     }
        !           864:     else
        !           865:     {
        !           866:       if (DEBUGLEVEL) fprintferr("candidate = %Z\n",CSF);
        !           867:       ESF = embedding_of_potential_subfields(CSF,DATA,listdelta);
        !           868:       if (!ESF) { avma = av2; continue; }
        !           869:
        !           870:       if (DEBUGLEVEL) fprintferr("embedding = %Z\n",ESF);
        !           871:       LSB = gerepileupto(av2, concat(LSB, _subfield(CSF,ESF)));
        !           872:     }
        !           873:   }
        !           874:   delete_var();
        !           875:   for (i=1; i<llist; i++) free((void*)listpotbl[i]);
        !           876:   free((void*)(listpotbl-1));
        !           877:   if (DEBUGLEVEL) fprintferr("\nSubfields of degree %ld: %Z\n",d, LSB);
        !           878:   return gerepilecopy(av, LSB);
        !           879: }
        !           880:
        !           881: static GEN
        !           882: fix_var(GEN x, long v)
        !           883: {
        !           884:   long i, l = lg(x);
        !           885:   x = gcopy(x); if (!v) return x;
        !           886:   for (i=1; i<l; i++) { GEN t = (GEN)x[i]; setvarn(t[1],v); setvarn(t[2],v); }
        !           887:   return x;
        !           888: }
        !           889:
        !           890: extern GEN get_nfpol(GEN x, GEN *nf);
        !           891: extern GEN vandermondeinverseprep(GEN T, GEN L);
        !           892: extern GEN ZX_disc_all(GEN x, ulong bound);
        !           893: extern GEN indexpartial(GEN P, GEN DP);
        !           894:
        !           895: void
        !           896: subfields_poldata(GEN T, struct poldata *PD)
        !           897: {
        !           898:   int     i, n;
        !           899:   GEN     L, z, prep, den;
        !           900:   long    prec;
        !           901:
        !           902:   GEN nf, dis;
        !           903:   T = get_nfpol(T, &nf);
        !           904:   PD->pol = dummycopy(T); /* may change variable number later */
        !           905:   if (nf)
        !           906:   {
        !           907:     PD->den = (GEN)nf[4];
        !           908:     PD->roo = (GEN)nf[6];
        !           909:     PD->dis = mulii(absi((GEN)nf[3]), sqri((GEN)nf[4]));
        !           910:     return;
        !           911:   }
        !           912:
        !           913:   /* copy-paste from galconj.c:galoisborne. Merge as soon as possible */
        !           914:   /* start of copy-paste */
        !           915:   n = degpol(T);
        !           916:   prec = 1;
        !           917:   for (i = 2; i < lgef(T); i++)
        !           918:     if (lg(T[i]) > prec)
        !           919:       prec = lg(T[i]);
        !           920:   /*prec++;*/
        !           921:   if (DEBUGLEVEL>=4) gentimer(3);
        !           922:   L = roots(T, prec);
        !           923:   if (DEBUGLEVEL>=4) genmsgtimer(3,"roots");
        !           924:   for (i = 1; i <= n; i++)
        !           925:   {
        !           926:     z = (GEN) L[i];
        !           927:     if (signe(z[2]))
        !           928:       break;
        !           929:     L[i] = z[1];
        !           930:   }
        !           931:   if (DEBUGLEVEL>=4) gentimer(3);
        !           932:   prep=vandermondeinverseprep(T, L);
        !           933:   /* end of copy-paste */
        !           934:   {
        !           935:     GEN res = divide_conquer_prod(gabs(prep,prec), mpmul);
        !           936:     disable_dbg(0);
        !           937:     dis = ZX_disc_all(T, 1+logint(res,gdeux,NULL));
        !           938:     disable_dbg(-1);
        !           939:     den = indexpartial(T,dis);
        !           940:   }
        !           941:
        !           942:   PD->den = den;
        !           943:   PD->roo = L;
        !           944:   PD->dis = absi(dis);
        !           945: }
        !           946:
        !           947: GEN
        !           948: subfields(GEN nf,GEN d)
        !           949: {
        !           950:   ulong av = avma;
        !           951:   long di,N,v0;
        !           952:   GEN LSB,pol;
        !           953:   struct poldata PD;
        !           954:
        !           955:   pol = get_nfpol(nf, &nf); /* in order to treat trivial cases */
        !           956:   v0=varn(pol); N=degpol(pol); di=itos(d);
        !           957:   if (di==N) return gerepilecopy(av, _subfield(pol, polx[v0]));
        !           958:   if (di==1) return gerepilecopy(av, _subfield(polx[v0], pol));
        !           959:   if (di < 1 || di > N || N % di) return cgetg(1,t_VEC);
        !           960:
        !           961:   subfields_poldata(nf, &PD);
        !           962:   pol = PD.pol;
        !           963:   setvarn(pol, 0);
        !           964:   LSB = subfields_of_given_degree(PD, di);
        !           965:   return gerepileupto(av, fix_var(LSB,v0));
        !           966: }
        !           967:
        !           968: static GEN
        !           969: subfieldsall(GEN nf)
        !           970: {
        !           971:   ulong av = avma;
        !           972:   long N,ld,i,v0;
        !           973:   GEN pol,dg,LSB,NLSB;
        !           974:   struct poldata PD;
        !           975:
        !           976:   subfields_poldata(nf, &PD);
        !           977:   pol = PD.pol;
        !           978:
        !           979:   v0 = varn(pol); N = degpol(pol);
        !           980:   dg = divisors(stoi(N)); ld = lg(dg)-1;
        !           981:   if (DEBUGLEVEL) fprintferr("\n***** Entering subfields\n\npol = %Z\n",pol);
        !           982:
        !           983:   LSB = _subfield(pol,polx[0]);
        !           984:   if (ld > 2)
        !           985:   {
        !           986:     setvarn(pol, 0);
        !           987:     for (i=2; i<ld; i++)
        !           988:     {
        !           989:       ulong av1 = avma;
        !           990:       NLSB = subfields_of_given_degree(PD, N / itos((GEN)dg[i]));
        !           991:       if (lg(NLSB) > 1) LSB = concatsp(LSB,NLSB); else avma = av1;
        !           992:     }
        !           993:   }
        !           994:   LSB = concatsp(LSB, _subfield(polx[0],pol));
        !           995:   if (DEBUGLEVEL) fprintferr("\n***** Leaving subfields\n\n");
        !           996:   return gerepileupto(av, fix_var(LSB,v0));
        !           997: }
        !           998:
        !           999: GEN
        !          1000: subfields0(GEN nf,GEN d)
        !          1001: {
        !          1002:   return d? subfields(nf,d): subfieldsall(nf);
        !          1003: }
        !          1004:
        !          1005: /* irreducible (unitary) polynomial of degree n over Fp[v] */
        !          1006: GEN
        !          1007: ffinit(GEN p,long n,long v)
        !          1008: {
        !          1009:   long av = avma;
        !          1010:   GEN pol;
        !          1011:
        !          1012:   if (n<=0) err(talker,"non positive degree in ffinit");
        !          1013:   if (typ(p) != t_INT) err(typeer,"ffinit");
        !          1014:   if (v<0) v = 0;
        !          1015:   for(;; avma = av)
        !          1016:   {
        !          1017:     pol = gadd(gpowgs(polx[v],n), FpX_rand(n-1,v, p));
        !          1018:     if (FpX_is_irred(pol, p)) break;
        !          1019:   }
        !          1020:   return gerepileupto(av, FpX(pol,p));
        !          1021: }

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