Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - modules - subfield.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24988-2584e74448) Lines: 913 941 97.0 %
Date: 2020-01-26 05:57:03 Functions: 47 47 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000-2004  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*               SUBFIELDS OF A NUMBER FIELD                       */
      17             : /*   J. Klueners and M. Pohst, J. Symb. Comp. (1996), vol. 11      */
      18             : /*                                                                 */
      19             : /*******************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : typedef struct _poldata {
      24             :   GEN pol;
      25             :   GEN dis; /* |disc(pol)| */
      26             :   GEN roo; /* roots(pol) */
      27             :   GEN den; /* multiple of index(pol) */
      28             : } poldata;
      29             : typedef struct _primedata {
      30             :   GEN p;  /* prime */
      31             :   GEN pol; /* pol mod p, squarefree */
      32             :   GEN ff; /* factorization of pol mod p */
      33             :   GEN Z; /* cycle structure of the above [ Frobenius orbits ] */
      34             :   long lcm; /* lcm of the above */
      35             :   GEN T;  /* ffinit(p, lcm) */
      36             : 
      37             :   GEN fk;      /* factorization of pol over F_(p^lcm) */
      38             :   GEN firstroot; /* *[i] = index of first root of fk[i] */
      39             :   GEN interp;    /* *[i] = interpolation polynomial for fk[i]
      40             :                   * [= 1 on the first root firstroot[i], 0 on the others] */
      41             :   GEN bezoutC; /* Bezout coefficients attached to the ff[i] */
      42             :   GEN Trk;     /* used to compute traces (cf poltrace) */
      43             : } primedata;
      44             : typedef struct _blockdata {
      45             :   poldata *PD; /* data depending from pol */
      46             :   primedata *S;/* data depending from pol, p */
      47             :   GEN DATA; /* data depending from pol, p, degree, # translations [updated] */
      48             :   long N; /* deg(PD.pol) */
      49             :   long d; /* subfield degree */
      50             :   long size;/* block degree = N/d */
      51             : } blockdata;
      52             : 
      53             : static GEN print_block_system(blockdata *B, GEN Y, GEN BS);
      54             : static GEN test_block(blockdata *B, GEN L, GEN D);
      55             : 
      56             : /* COMBINATORIAL PART: generate potential block systems */
      57             : 
      58             : #define BIL 32 /* for 64bit machines also */
      59             : /* Computation of potential block systems of given size d attached to a
      60             :  * rational prime p: give a row vector of row vectors containing the
      61             :  * potential block systems of imprimitivity; a potential block system is a
      62             :  * vector of row vectors (enumeration of the roots). */
      63             : static GEN
      64         630 : calc_block(blockdata *B, GEN Z, GEN Y, GEN SB)
      65             : {
      66         630 :   long r = lg(Z), lK, i, j, t, tp, T, u, nn, lnon, lY;
      67             :   GEN K, n, non, pn, pnon, e, Yp, Zp, Zpp;
      68         630 :   pari_sp av0 = avma;
      69             : 
      70         630 :   if (DEBUGLEVEL>3)
      71             :   {
      72           0 :     err_printf("lg(Z) = %ld, lg(Y) = %ld\n", r,lg(Y));
      73           0 :     if (DEBUGLEVEL > 5)
      74             :     {
      75           0 :       err_printf("Z = %Ps\n",Z);
      76           0 :       err_printf("Y = %Ps\n",Y);
      77             :     }
      78             :   }
      79         630 :   lnon = minss(BIL, r);
      80         630 :   e    = new_chunk(BIL);
      81         630 :   n    = new_chunk(r);
      82         630 :   non  = new_chunk(lnon);
      83         630 :   pnon = new_chunk(lnon);
      84         630 :   pn   = new_chunk(lnon);
      85             : 
      86         630 :   Zp   = cgetg(lnon,t_VEC);
      87         630 :   Zpp  = cgetg(lnon,t_VEC); nn = 0;
      88         630 :   for (i=1; i<r; i++) { n[i] = lg(gel(Z,i))-1; nn += n[i]; }
      89         630 :   lY = lg(Y); Yp = cgetg(lY+1,t_VEC);
      90         630 :   for (j=1; j<lY; j++) gel(Yp,j) = gel(Y,j);
      91             : 
      92             :   {
      93         630 :     pari_sp av = avma;
      94         630 :     long k = nn / B->size;
      95        1309 :     for (j = 1; j < r; j++)
      96         777 :       if (n[j] % k) break;
      97         630 :     if (j == r)
      98             :     {
      99         532 :       gel(Yp,lY) = Z;
     100         532 :       SB = print_block_system(B, Yp, SB);
     101         532 :       set_avma(av);
     102             :     }
     103             :   }
     104         630 :   gel(Yp,lY) = Zp;
     105             : 
     106         630 :   K = divisorsu(n[1]); lK = lg(K);
     107        3087 :   for (i=1; i<lK; i++)
     108             :   {
     109        2457 :     long ngcd = n[1], k = K[i], dk = B->size*k, lpn = 0;
     110        3381 :     for (j=2; j<r; j++)
     111         924 :       if (n[j]%k == 0)
     112             :       {
     113         868 :         if (++lpn >= BIL) pari_err_OVERFLOW("calc_block");
     114         868 :         pn[lpn] = n[j]; pnon[lpn] = j;
     115         868 :         ngcd = ugcd(ngcd, n[j]);
     116             :       }
     117        2457 :     if (dk % ngcd) continue;
     118        1505 :     T = 1L<<lpn;
     119        1505 :     if (lpn == r-2)
     120             :     {
     121        1449 :       T--; /* done already above --> print_block_system */
     122        1449 :       if (!T) continue;
     123             :     }
     124             : 
     125         525 :     if (dk == n[1])
     126             :     { /* empty subset, t = 0. Split out for clarity */
     127         140 :       Zp[1] = Z[1]; setlg(Zp, 2);
     128         140 :       for (u=1,j=2; j<r; j++) Zpp[u++] = Z[j];
     129         140 :       setlg(Zpp, u);
     130         140 :       SB = calc_block(B, Zpp, Yp, SB);
     131             :     }
     132             : 
     133         721 :     for (t = 1; t < T; t++)
     134             :     { /* loop through all non-empty subsets of [1..lpn] */
     135         714 :       for (nn=n[1],tp=t, u=1; u<=lpn; u++,tp>>=1)
     136             :       {
     137         518 :         if (tp&1) { nn += pn[u]; e[u] = 1; } else e[u] = 0;
     138             :       }
     139         196 :       if (dk != nn) continue;
     140             : 
     141          35 :       for (j=1; j<r; j++) non[j]=0;
     142          35 :       Zp[1] = Z[1];
     143         126 :       for (u=2,j=1; j<=lpn; j++)
     144          91 :         if (e[j]) { Zp[u] = Z[pnon[j]]; non[pnon[j]] = 1; u++; }
     145          35 :       setlg(Zp, u);
     146         126 :       for (u=1,j=2; j<r; j++)
     147          91 :         if (!non[j]) Zpp[u++] = Z[j];
     148          35 :       setlg(Zpp, u);
     149          35 :       SB = calc_block(B, Zpp, Yp, SB);
     150             :     }
     151             :   }
     152         630 :   return gc_const(av0, SB);
     153             : }
     154             : 
     155             : /* product of permutations. Put the result in perm1. */
     156             : static void
     157      208061 : perm_mul_i(GEN perm1, GEN perm2)
     158             : {
     159      208061 :   long i, N = lg(perm1);
     160      208061 :   pari_sp av = avma;
     161      208061 :   GEN perm = new_chunk(N);
     162      208061 :   for (i=1; i<N; i++) perm[i] = perm1[perm2[i]];
     163      208061 :   for (i=1; i<N; i++) perm1[i]= perm[i];
     164      208061 :   set_avma(av);
     165      208061 : }
     166             : 
     167             : /* cy is a cycle; compute cy^l as a permutation */
     168             : static GEN
     169       38486 : cycle_power_to_perm(GEN perm,GEN cy,long l)
     170             : {
     171       38486 :   long lp,i,j,b, N = lg(perm), lcy = lg(cy)-1;
     172             : 
     173       38486 :   lp = l % lcy;
     174       38486 :   for (i=1; i<N; i++) perm[i] = i;
     175       38486 :   if (lp)
     176             :   {
     177       35042 :     pari_sp av = avma;
     178       35042 :     GEN p1 = new_chunk(N);
     179       35042 :     b = cy[1];
     180       35042 :     for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
     181       35042 :     perm[b] = cy[1];
     182       35042 :     for (i=1; i<N; i++) p1[i] = perm[i];
     183             : 
     184       35042 :     for (j=2; j<=lp; j++) perm_mul_i(perm,p1);
     185       35042 :     set_avma(av);
     186             :   }
     187       38486 :   return perm;
     188             : }
     189             : 
     190             : /* image du block system D par la permutation perm */
     191             : static GEN
     192       13475 : im_block_by_perm(GEN D,GEN perm)
     193             : {
     194       13475 :   long i, lb = lg(D);
     195       13475 :   GEN Dn = cgetg(lb,t_VEC);
     196       13475 :   for (i=1; i<lb; i++) gel(Dn,i) = vecsmallpermute(perm, gel(D,i));
     197       13475 :   return Dn;
     198             : }
     199             : 
     200             : static void
     201         756 : append(GEN D, GEN a)
     202             : {
     203         756 :   long i,l = lg(D), m = lg(a);
     204         756 :   GEN x = D + (l-1);
     205         756 :   for (i=1; i<m; i++) gel(x,i) = gel(a,i);
     206         756 :   setlg(D, l+m-1);
     207         756 : }
     208             : 
     209             : static GEN
     210         532 : print_block_system(blockdata *B, GEN Y, GEN SB)
     211             : {
     212         532 :   long i, j, l, ll, lp, u, v, ns, r = lg(Y), N = B->N;
     213             :   long *k, *n, **e, *t;
     214         532 :   GEN D, De, Z, cyperm, perm, VOID = cgetg(1, t_VECSMALL);
     215             : 
     216         532 :   if (DEBUGLEVEL>5) err_printf("Y = %Ps\n",Y);
     217         532 :   n = new_chunk(N+1);
     218         532 :   D = vectrunc_init(N+1);
     219         532 :   t = new_chunk(r+1);
     220         532 :   k = new_chunk(r+1);
     221         532 :   Z = cgetg(r+1, t_VEC);
     222        1288 :   for (ns=0,i=1; i<r; i++)
     223             :   {
     224         756 :     GEN Yi = gel(Y,i);
     225         756 :     long ki = 0, si = lg(Yi)-1;
     226             : 
     227         756 :     for (j=1; j<=si; j++) { n[j] = lg(gel(Yi,j))-1; ki += n[j]; }
     228         756 :     ki /= B->size;
     229         756 :     De = cgetg(ki+1,t_VEC);
     230         756 :     for (j=1; j<=ki; j++) gel(De,j) = VOID;
     231        1715 :     for (j=1; j<=si; j++)
     232             :     {
     233         959 :       GEN cy = gel(Yi,j);
     234        9303 :       for (l=1,lp=0; l<=n[j]; l++)
     235             :       {
     236        8344 :         lp++; if (lp > ki) lp = 1;
     237        8344 :         gel(De,lp) = vecsmall_append(gel(De,lp), cy[l]);
     238             :       }
     239             :     }
     240         756 :     append(D, De);
     241         756 :     if (si>1 && ki>1)
     242             :     {
     243         182 :       GEN p1 = cgetg(si,t_VEC);
     244         182 :       for (j=2; j<=si; j++) p1[j-1] = Yi[j];
     245         182 :       ns++;
     246         182 :       t[ns] = si-1;
     247         182 :       k[ns] = ki-1;
     248         182 :       gel(Z,ns) = p1;
     249             :     }
     250             :   }
     251         532 :   if (DEBUGLEVEL>2) err_printf("\nns = %ld\n",ns);
     252         532 :   if (!ns) return test_block(B, SB, D);
     253             : 
     254         161 :   setlg(Z, ns+1);
     255         161 :   e = (long**)new_chunk(ns+1);
     256         343 :   for (i=1; i<=ns; i++)
     257             :   {
     258         182 :     e[i] = new_chunk(t[i]+1);
     259         182 :     for (j=1; j<=t[i]; j++) e[i][j] = 0;
     260             :   }
     261         161 :   cyperm= cgetg(N+1,t_VECSMALL);
     262         161 :   perm  = cgetg(N+1,t_VECSMALL); i = ns;
     263             :   do
     264             :   {
     265       13475 :     pari_sp av = avma;
     266       13475 :     for (u=1; u<=N; u++) perm[u] = u;
     267       27706 :     for (u=1; u<=ns; u++)
     268       52717 :       for (v=1; v<=t[u]; v++)
     269       38486 :         perm_mul_i(perm, cycle_power_to_perm(cyperm, gmael(Z,u,v), e[u][v]));
     270       13475 :     SB = test_block(B, SB, im_block_by_perm(D,perm));
     271       13475 :     set_avma(av);
     272             : 
     273             :     /* i = 1..ns, j = 1..t[i], e[i][j] loop through 0..k[i].
     274             :      * TODO: flatten to 1-dimensional loop */
     275       13475 :     if (++e[ns][t[ns]] > k[ns])
     276             :     {
     277        1281 :       j = t[ns]-1;
     278        1281 :       while (j>=1 && e[ns][j] == k[ns]) j--;
     279        1281 :       if (j >= 1) { e[ns][j]++; for (l=j+1; l<=t[ns]; l++) e[ns][l] = 0; }
     280             :       else
     281             :       {
     282         266 :         i = ns-1;
     283         553 :         while (i>=1)
     284             :         {
     285         126 :           j = t[i];
     286         126 :           while (j>=1 && e[i][j] == k[i]) j--;
     287         126 :           if (j<1) i--;
     288             :           else
     289             :           {
     290         105 :             e[i][j]++;
     291         105 :             for (l=j+1; l<=t[i]; l++) e[i][l] = 0;
     292         210 :             for (ll=i+1; ll<=ns; ll++)
     293         105 :               for (l=1; l<=t[ll]; l++) e[ll][l] = 0;
     294         105 :             break;
     295             :           }
     296             :         }
     297             :       }
     298             :     }
     299             :   }
     300       13475 :   while (i > 0);
     301         161 :   return SB;
     302             : }
     303             : 
     304             : /* ALGEBRAIC PART: test potential block systems */
     305             : 
     306             : static GEN
     307        5110 : polsimplify(GEN x)
     308             : {
     309        5110 :   long i,lx = lg(x);
     310       67081 :   for (i=2; i<lx; i++)
     311       61971 :     if (typ(gel(x,i)) == t_POL) gel(x,i) = constant_coeff(gel(x,i));
     312        5110 :   return x;
     313             : }
     314             : 
     315             : /* return 0 if |g[i]| > M[i] for some i; 1 otherwise */
     316             : static long
     317        5110 : ok_coeffs(GEN g,GEN M)
     318             : {
     319        5110 :   long i, lg = lg(g)-1; /* g is monic, and cst term is ok */
     320        7441 :   for (i=3; i<lg; i++)
     321        6825 :     if (abscmpii(gel(g,i), gel(M,i)) > 0) return 0;
     322         616 :   return 1;
     323             : }
     324             : 
     325             : /* assume x in Fq, return Tr_{Fq/Fp}(x) as a t_INT */
     326             : static GEN
     327        9107 : trace(GEN x, GEN Trq, GEN p)
     328             : {
     329             :   long i, l;
     330             :   GEN s;
     331        9107 :   if (typ(x) == t_INT) return Fp_mul(x, gel(Trq,1), p);
     332        9107 :   l = lg(x)-1; if (l == 1) return gen_0;
     333        9107 :   x++; s = mulii(gel(x,1), gel(Trq,1));
     334       94605 :   for (i=2; i<l; i++)
     335       85498 :     s = addii(s, mulii(gel(x,i), gel(Trq,i)));
     336        9107 :   return modii(s, p);
     337             : }
     338             : 
     339             : /* assume x in Fq[X], return Tr_{Fq[X]/Fp[X]}(x), varn(X) = 0 */
     340             : static GEN
     341        1050 : poltrace(GEN x, GEN Trq, GEN p)
     342             : {
     343             :   long i,l;
     344             :   GEN y;
     345        1050 :   if (typ(x) == t_INT || varn(x) != 0) return trace(x, Trq, p);
     346        1050 :   y = cgetg_copy(x, &l); y[1] = x[1];
     347        1050 :   for (i=2; i<l; i++) gel(y,i) = trace(gel(x,i),Trq,p);
     348        1050 :   return normalizepol_lg(y, l);
     349             : }
     350             : 
     351             : /* Find h in Fp[X] such that h(a[i]) = listdelta[i] for all modular factors
     352             :  * ff[i], where a[i] is a fixed root of ff[i] in Fq = Z[Y]/(p,T) [namely the
     353             :  * first one in FpX_factorff_irred output]. Let f = ff[i], A the given root,
     354             :  * then h mod f is Tr_Fq/Fp ( h(A) f(X)/(X-A)f'(A) ), most of the expression
     355             :  * being precomputed. The complete h is recovered via chinese remaindering */
     356             : static GEN
     357         532 : chinese_retrieve_pol(GEN DATA, primedata *S, GEN listdelta)
     358             : {
     359         532 :   GEN interp, bezoutC, h, p = S->p, pol = FpX_red(gel(DATA,1), p);
     360             :   long i, l;
     361         532 :   interp = gel(DATA,9);
     362         532 :   bezoutC= gel(DATA,6);
     363             : 
     364         532 :   h = NULL; l = lg(interp);
     365        1582 :   for (i=1; i<l; i++)
     366             :   { /* h(firstroot[i]) = listdelta[i] */
     367        1050 :     GEN t = FqX_Fq_mul(gel(interp,i), gel(listdelta,i), S->T, p);
     368        1050 :     t = poltrace(t, gel(S->Trk,i), p);
     369        1050 :     t = FpX_mul(t, gel(bezoutC,i), p);
     370        1050 :     h = h? FpX_add(h,t,p): t;
     371             :   }
     372         532 :   return FpX_rem(h, pol, p);
     373             : }
     374             : 
     375             : /* g in Z[X] potentially defines a subfield of Q[X]/f. It is a subfield iff A
     376             :  * (cf subfield) was a block system; then there
     377             :  * exists h in Q[X] such that f | g o h. listdelta determines h s.t f | g o h
     378             :  * in Fp[X] (cf chinese_retrieve_pol). Try to lift it; den is a
     379             :  * multiplicative bound for denominator of lift. */
     380             : static GEN
     381         532 : embedding(GEN g, GEN DATA, primedata *S, GEN den, GEN listdelta)
     382             : {
     383         532 :   GEN TR, w0_Q, w0, w1_Q, w1, wpow, h0, gp, T, q2, q, maxp, a, p = S->p;
     384             :   long rt;
     385             :   pari_sp av;
     386             : 
     387         532 :   T   = gel(DATA,1); rt = brent_kung_optpow(degpol(T), 4, 3);
     388         532 :   maxp= gel(DATA,7);
     389         532 :   gp = RgX_deriv(g); av = avma;
     390         532 :   w0 = chinese_retrieve_pol(DATA, S, listdelta);
     391         532 :   w0_Q = centermod(gmul(w0,den), p);
     392         532 :   h0 = FpXQ_inv(FpX_FpXQ_eval(gp,w0, T,p), T,p); /* = 1/g'(w0) mod (T,p) */
     393         532 :   wpow = NULL; q = sqri(p);
     394             :   for(;;)
     395             :   {/* Given g,w0,h0 in Z[x], s.t. h0.g'(w0) = 1 and g(w0) = 0 mod (T,p), find
     396             :     * [w1,h1] satisfying the same conditions mod p^2, [w1,h1] = [w0,h0] (mod p)
     397             :     * (cf. Dixon: J. Austral. Math. Soc., Series A, vol.49, 1990, p.445) */
     398        3668 :     if (DEBUGLEVEL>1)
     399           0 :       err_printf("lifting embedding mod p^k = %Ps^%ld\n",S->p, Z_pval(q,S->p));
     400             : 
     401             :     /* w1 := w0 - h0 g(w0) mod (T,q) */
     402        2100 :     if (wpow) a = FpX_FpXQV_eval(g,wpow, T,q);
     403         532 :     else      a = FpX_FpXQ_eval(g,w0, T,q); /* first time */
     404             :     /* now, a = 0 (p) */
     405        2100 :     a = FpXQ_mul(ZX_neg(h0), ZX_Z_divexact(a, p), T,p);
     406        2100 :     w1 = ZX_add(w0, ZX_Z_mul(a, p));
     407             : 
     408        2100 :     w1_Q = centermod(ZX_Z_mul(w1, remii(den,q)), q);
     409        2100 :     if (ZX_equal(w1_Q, w0_Q))
     410             :     {
     411         392 :       GEN G = is_pm1(den)? g: RgX_rescale(g,den);
     412         392 :       if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
     413             :     }
     414        1708 :     else if (cmpii(q,maxp) > 0)
     415             :     {
     416         140 :       GEN G = is_pm1(den)? g: RgX_rescale(g,den);
     417         140 :       if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
     418           0 :       if (DEBUGLEVEL) err_printf("coeff too big for embedding\n");
     419           0 :       return NULL;
     420             :     }
     421        1568 :     gerepileall(av, 5, &w1,&h0,&w1_Q,&q,&p);
     422        1568 :     q2 = sqri(q);
     423        1568 :     wpow = FpXQ_powers(w1, rt, T, q2);
     424             :     /* h0 := h0 * (2 - h0 g'(w1)) mod (T,q)
     425             :      *     = h0 + h0 * (1 - h0 g'(w1)) */
     426        1568 :     a = FpXQ_mul(ZX_neg(h0), FpX_FpXQV_eval(gp, FpXV_red(wpow,q),T,q), T,q);
     427        1568 :     a = ZX_Z_add_shallow(a, gen_1); /* 1 - h0 g'(w1) = 0 (p) */
     428        1568 :     a = FpXQ_mul(h0, ZX_Z_divexact(a, p), T,p);
     429        1568 :     h0 = ZX_add(h0, ZX_Z_mul(a, p));
     430        1568 :     w0 = w1; w0_Q = w1_Q; p = q; q = q2;
     431             :   }
     432         532 :   TR = gel(DATA,5);
     433         532 :   if (!gequal0(TR)) w1_Q = RgX_translate(w1_Q, TR);
     434         532 :   return gdiv(w1_Q,den);
     435             : }
     436             : 
     437             : /* return U list of polynomials s.t U[i] = 1 mod fk[i] and 0 mod fk[j] for all
     438             :  * other j */
     439             : static GEN
     440         266 : get_bezout(GEN pol, GEN fk, GEN p)
     441             : {
     442         266 :   long i, l = lg(fk);
     443         266 :   GEN A, B, d, u, v, U = cgetg(l, t_VEC);
     444         623 :   for (i=1; i<l; i++)
     445             :   {
     446         357 :     A = gel(fk,i);
     447         357 :     B = FpX_div(pol, A, p);
     448         357 :     d = FpX_extgcd(A,B,p, &u, &v);
     449         357 :     if (degpol(d) > 0) pari_err_COPRIME("get_bezout",A,B);
     450         357 :     d = constant_coeff(d);
     451         357 :     if (!gequal1(d)) v = FpX_Fp_mul(v, Fp_inv(d, p), p);
     452         357 :     gel(U,i) = FpX_mul(B,v, p);
     453             :   }
     454         266 :   return U;
     455             : }
     456             : 
     457             : static GEN
     458         266 : init_traces(GEN ff, GEN T, GEN p)
     459             : {
     460         266 :   long N = degpol(T),i,j,k, r = lg(ff);
     461         266 :   GEN Frob = FpX_matFrobenius(T,p);
     462             :   GEN y,p1,p2,Trk,pow,pow1;
     463             : 
     464         266 :   k = degpol(gel(ff,r-1)); /* largest degree in modular factorization */
     465         266 :   pow = cgetg(k+1, t_VEC);
     466         266 :   gel(pow,1) = gen_0; /* dummy */
     467         266 :   gel(pow,2) = Frob;
     468         266 :   pow1= cgetg(k+1, t_VEC); /* 1st line */
     469        1575 :   for (i=3; i<=k; i++)
     470        1309 :     gel(pow,i) = FpM_mul(gel(pow,i-1), Frob, p);
     471         266 :   gel(pow1,1) = gen_0; /* dummy */
     472        1841 :   for (i=2; i<=k; i++)
     473             :   {
     474        1575 :     p1 = cgetg(N+1, t_VEC);
     475        1575 :     gel(pow1,i) = p1; p2 = gel(pow,i);
     476        1575 :     for (j=1; j<=N; j++) gel(p1,j) = gcoeff(p2,1,j);
     477             :   }
     478             : 
     479             :   /* Trk[i] = line 1 of x -> x + x^p + ... + x^{p^(i-1)} */
     480         266 :   Trk = pow; /* re-use (destroy) pow */
     481         266 :   gel(Trk,1) = vec_ei(N,1);
     482        1841 :   for (i=2; i<=k; i++)
     483        1575 :     gel(Trk,i) = gadd(gel(Trk,i-1), gel(pow1,i));
     484         266 :   y = cgetg(r, t_VEC);
     485         266 :   for (i=1; i<r; i++) y[i] = Trk[degpol(gel(ff,i))];
     486         266 :   return y;
     487             : }
     488             : 
     489             : static void
     490         266 : init_primedata(primedata *S)
     491             : {
     492         266 :   long i, j, l, lff = lg(S->ff), v = fetch_var(), N = degpol(S->pol);
     493         266 :   GEN T, p = S->p;
     494             : 
     495         266 :   if (S->lcm == degpol(gel(S->ff,lff-1)))
     496             :   {
     497         252 :     T = leafcopy(gel(S->ff,lff-1));
     498         252 :     setvarn(T, v);
     499             :   }
     500             :   else
     501          14 :     T = init_Fq(p, S->lcm, v);
     502         266 :   S->T = T;
     503         266 :   S->firstroot = cgetg(lff, t_VECSMALL);
     504         266 :   S->interp = cgetg(lff, t_VEC);
     505         266 :   S->fk = cgetg(N+1, t_VEC);
     506         623 :   for (l=1,j=1; j<lff; j++)
     507             :   { /* compute roots and fix ordering (Frobenius cycles) */
     508         357 :     GEN F = gel(S->ff, j), deg1 = FpX_factorff_irred(F, T,p);
     509         357 :     GEN H = gel(deg1,1), a = Fq_neg(constant_coeff(H), T,p);
     510         357 :     GEN Q = FqX_div(F, H, T,p);
     511         357 :     GEN q = Fq_inv(FqX_eval(Q, a, T,p), T,p);
     512         357 :     gel(S->interp,j) = FqX_Fq_mul(Q, q, T,p); /* = 1 at a, 0 at other roots */
     513         357 :     S->firstroot[j] = l;
     514         357 :     for (i=1; i<lg(deg1); i++,l++) gel(S->fk, l) = gel(deg1, i);
     515             :   }
     516         266 :   S->Trk     = init_traces(S->ff, T,p);
     517         266 :   S->bezoutC = get_bezout(S->pol, S->ff, p);
     518         266 : }
     519             : 
     520             : static int
     521         273 : choose_prime(primedata *S, GEN pol)
     522             : {
     523         273 :   long i, j, k, r, lcm, oldr, K, N = degpol(pol);
     524             :   ulong p, pp;
     525             :   GEN Z, ff, n, oldn;
     526             :   pari_sp av;
     527             :   forprime_t T;
     528             : 
     529         273 :   u_forprime_init(&T, (N*N) >> 2, ULONG_MAX);
     530         273 :   oldr = S->lcm = LONG_MAX;
     531         273 :   S->ff = oldn = NULL; pp = 0; /* gcc -Wall */
     532         273 :   av = avma; K = N + 10;
     533        2548 :   for(k = 1; k < K || !S->ff; k++,set_avma(av))
     534             :   {
     535             :     GEN Tp;
     536        2471 :     if (k > 5 * N) return 0;
     537             :     do
     538             :     {
     539        2562 :       p = u_forprime_next(&T);
     540        2562 :       Tp = ZX_to_Flx(pol, p);
     541             :     }
     542        2562 :     while (!Flx_is_squarefree(Tp, p));
     543        2471 :     ff = gel(Flx_factor(Tp, p), 1);
     544        2471 :     r = lg(ff)-1;
     545        2471 :     if (r == N || r >= BIL) continue;
     546             : 
     547        2436 :     n = cgetg(r+1, t_VECSMALL); lcm = n[1] = degpol(gel(ff,1));
     548        2436 :     for (j=2; j<=r; j++) { n[j] = degpol(gel(ff,j)); lcm = ulcm(lcm, n[j]); }
     549        2436 :     if (r > oldr || (r == oldr && (lcm <= S->lcm || S->lcm > 2*N)))
     550        2002 :       continue;
     551         434 :     if (DEBUGLEVEL) err_printf("p = %lu,\tlcm = %ld,\torbits: %Ps\n",p,lcm,n);
     552             : 
     553         434 :     pp = p;
     554         434 :     oldr = r;
     555         434 :     oldn = n;
     556         434 :     S->ff = ff;
     557         434 :     S->lcm = lcm; if (r == 1) break;
     558         238 :     av = avma;
     559             :   }
     560         273 :   if (oldr > 6) return 0;
     561         266 :   if (DEBUGLEVEL) err_printf("Chosen prime: p = %ld\n", pp);
     562         266 :   FlxV_to_ZXV_inplace(S->ff);
     563         266 :   S->p  = utoipos(pp);
     564         266 :   S->pol = FpX_red(pol, S->p); init_primedata(S);
     565         266 :   n = oldn; r = lg(n); S->Z = Z = cgetg(r,t_VEC);
     566         623 :   for (k=0,i=1; i<r; i++)
     567             :   {
     568         357 :     GEN t = cgetg(n[i]+1, t_VECSMALL); gel(Z,i) = t;
     569         357 :     for (j=1; j<=n[i]; j++) t[j] = ++k;
     570             :   }
     571         266 :   return 1;
     572             : }
     573             : 
     574             : /* maxroot t_REAL */
     575             : static GEN
     576         539 : bound_for_coeff(long m, GEN rr, GEN *maxroot)
     577             : {
     578         539 :   long i,r1, lrr=lg(rr);
     579         539 :   GEN p1,b1,b2,B,M, C = matpascal(m-1);
     580             : 
     581        1687 :   for (r1=1; r1 < lrr; r1++)
     582        1673 :     if (typ(gel(rr,r1)) != t_REAL) break;
     583         539 :   r1--;
     584             : 
     585         539 :   rr = gabs(rr,0); *maxroot = vecmax(rr);
     586        6748 :   for (i=1; i<lrr; i++)
     587        6209 :     if (gcmp(gel(rr,i), gen_1) < 0) gel(rr,i) = gen_1;
     588         539 :   for (b1=gen_1,i=1; i<=r1; i++) b1 = gmul(b1, gel(rr,i));
     589         539 :   for (b2=gen_1    ; i<lrr; i++) b2 = gmul(b2, gel(rr,i));
     590         539 :   B = gmul(b1, gsqr(b2)); /* Mahler measure */
     591         539 :   M = cgetg(m+2, t_VEC); gel(M,1) = gel(M,2) = gen_0; /* unused */
     592        1988 :   for (i=1; i<m; i++)
     593             :   {
     594        1449 :     p1 = gadd(gmul(gcoeff(C, m, i+1), B),/* binom(m-1, i)   */
     595        1449 :               gcoeff(C, m, i));          /* binom(m-1, i-1) */
     596        1449 :     gel(M,i+2) = ceil_safe(p1);
     597             :   }
     598         539 :   return M;
     599             : }
     600             : 
     601             : /* d = requested degree for subfield. Return DATA, valid for given pol, S and d
     602             :  * If DATA != NULL, translate pol [ --> pol(X+1) ] and update DATA
     603             :  * 1: polynomial pol
     604             :  * 2: p^e (for Hensel lifts) such that p^e > max(M),
     605             :  * 3: Hensel lift to precision p^e of DATA[4]
     606             :  * 4: roots of pol in F_(p^S->lcm),
     607             :  * 5: number of polynomial changes (translations)
     608             :  * 6: Bezout coefficients attached to the S->ff[i]
     609             :  * 7: Hadamard bound for coefficients of h(x) such that g o h = 0 mod pol.
     610             :  * 8: bound M for polynomials defining subfields x PD->den
     611             :  * 9: *[i] = interpolation polynomial for S->ff[i] [= 1 on the first root
     612             :       S->firstroot[i], 0 on the others] */
     613             : static void
     614         539 : compute_data(blockdata *B)
     615             : {
     616             :   GEN ffL, roo, pe, p1, p2, fk, fhk, MM, maxroot, pol;
     617         539 :   primedata *S = B->S;
     618         539 :   GEN p = S->p, T = S->T, ff = S->ff, DATA = B->DATA;
     619         539 :   long i, j, l, e, N, lff = lg(ff);
     620             : 
     621         539 :   if (DEBUGLEVEL>1) err_printf("Entering compute_data()\n\n");
     622         539 :   pol = B->PD->pol; N = degpol(pol);
     623         539 :   roo = B->PD->roo;
     624         539 :   if (DATA)
     625             :   {
     626          84 :     GEN Xm1 = gsub(pol_x(varn(pol)), gen_1);
     627          84 :     GEN TR = addiu(gel(DATA,5), 1);
     628          84 :     GEN mTR = negi(TR), interp, bezoutC;
     629             : 
     630          84 :     if (DEBUGLEVEL>1) err_printf("... update (translate) an existing DATA\n\n");
     631             : 
     632          84 :     gel(DATA,5) = TR;
     633          84 :     pol = RgX_translate(gel(DATA,1), gen_m1);
     634          84 :     p1 = cgetg_copy(roo, &l);
     635          84 :     for (i=1; i<l; i++) gel(p1,i) = gadd(TR, gel(roo,i));
     636          84 :     roo = p1;
     637             : 
     638          84 :     fk = gel(DATA,4); l = lg(fk);
     639          84 :     for (i=1; i<l; i++) gel(fk,i) = gsub(Xm1, gel(fk,i));
     640             : 
     641          84 :     bezoutC = gel(DATA,6); l = lg(bezoutC);
     642          84 :     interp  = gel(DATA,9);
     643         210 :     for (i=1; i<l; i++)
     644             :     {
     645         126 :       if (degpol(gel(interp,i)) > 0) /* do not turn pol_1(0) into gen_1 */
     646             :       {
     647         126 :         p1 = RgX_translate(gel(interp,i), gen_m1);
     648         126 :         gel(interp,i) = FpXX_red(p1, p);
     649             :       }
     650         126 :       if (degpol(gel(bezoutC,i)) > 0)
     651             :       {
     652          77 :         p1 = RgX_translate(gel(bezoutC,i), gen_m1);
     653          77 :         gel(bezoutC,i) = FpXX_red(p1, p);
     654             :       }
     655             :     }
     656          84 :     ff = cgetg(lff, t_VEC); /* copy, do not overwrite! */
     657         210 :     for (i=1; i<lff; i++)
     658         126 :       gel(ff,i) = FpX_red(RgX_translate(gel(S->ff,i), mTR), p);
     659             :   }
     660             :   else
     661             :   {
     662         455 :     DATA = cgetg(10,t_VEC);
     663         455 :     fk = S->fk;
     664         455 :     gel(DATA,5) = gen_0;
     665         455 :     gel(DATA,6) = leafcopy(S->bezoutC);
     666         455 :     gel(DATA,9) = leafcopy(S->interp);
     667             :   }
     668         539 :   gel(DATA,1) = pol;
     669         539 :   MM = gmul2n(bound_for_coeff(B->d, roo, &maxroot), 1);
     670         539 :   gel(DATA,8) = MM;
     671         539 :   e = logintall(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [d-1 test] */
     672         539 :   gel(DATA,2) = pe;
     673         539 :   gel(DATA,4) = roots_from_deg1(fk);
     674             : 
     675             :   /* compute fhk = ZpX_liftfact(pol,fk,T,p,e,pe) in 2 steps
     676             :    * 1) lift in Zp to precision p^e */
     677         539 :   ffL = ZpX_liftfact(pol, ff, pe, p, e);
     678         539 :   fhk = NULL;
     679        1323 :   for (l=i=1; i<lff; i++)
     680             :   { /* 2) lift factorization of ff[i] in Qp[X] / T */
     681         784 :     GEN F, L = gel(ffL,i);
     682         784 :     long di = degpol(L);
     683         784 :     F = cgetg(di+1, t_VEC);
     684         784 :     for (j=1; j<=di; j++) F[j] = fk[l++];
     685         784 :     L = ZqX_liftfact(L, F, T, pe, p, e);
     686         784 :     fhk = fhk? shallowconcat(fhk, L): L;
     687             :   }
     688         539 :   gel(DATA,3) = roots_from_deg1(fhk);
     689             : 
     690         539 :   p1 = mulur(N, powruhalf(utor(N-1,DEFAULTPREC), N-1));
     691         539 :   p2 = powru(maxroot, B->size + N*(N-1)/2);
     692         539 :   p1 = divrr(mulrr(p1,p2), gsqrt(B->PD->dis,DEFAULTPREC));
     693         539 :   gel(DATA,7) = mulii(shifti(ceil_safe(p1), 1), B->PD->den);
     694             : 
     695         539 :   if (DEBUGLEVEL>1) {
     696           0 :     err_printf("f = %Ps\n",DATA[1]);
     697           0 :     err_printf("p = %Ps, lift to p^%ld\n", p, e);
     698           0 :     err_printf("2 * Hadamard bound * ind = %Ps\n",DATA[7]);
     699           0 :     err_printf("2 * M = %Ps\n",DATA[8]);
     700             :   }
     701         539 :   if (B->DATA) { DATA = gclone(DATA); clone_unlock(B->DATA); }
     702         539 :   B->DATA = DATA;
     703         539 : }
     704             : 
     705             : /* g = polynomial, h = embedding. Return [[g,h]] */
     706             : static GEN
     707        1463 : _subfield(GEN g, GEN h) { return mkvec(mkvec2(g,h)); }
     708             : 
     709             : /* Return a subfield, gen_0 [ change p ] or NULL [ not a subfield ] */
     710             : static GEN
     711       13930 : subfield(GEN A, blockdata *B)
     712             : {
     713       13930 :   long N, i, j, d, lf, m = lg(A)-1;
     714             :   GEN M, pe, pol, fhk, g, e, d_1_term, delta, listdelta, whichdelta;
     715       13930 :   GEN T = B->S->T, p = B->S->p, firstroot = B->S->firstroot;
     716             : 
     717       13930 :   pol= gel(B->DATA,1); N = degpol(pol); d = N/m; /* m | N */
     718       13930 :   pe = gel(B->DATA,2);
     719       13930 :   fhk= gel(B->DATA,3);
     720       13930 :   M  = gel(B->DATA,8);
     721             : 
     722       13930 :   delta = cgetg(m+1,t_VEC);
     723       13930 :   whichdelta = cgetg(N+1, t_VECSMALL);
     724       13930 :   d_1_term = gen_0;
     725      174202 :   for (i=1; i<=m; i++)
     726             :   {
     727      160272 :     GEN Ai = gel(A,i), p1 = gel(fhk,Ai[1]);
     728      637973 :     for (j=2; j<=d; j++)
     729      477701 :       p1 = Fq_mul(p1, gel(fhk,Ai[j]), T, pe);
     730      160272 :     gel(delta,i) = p1;
     731      160272 :     if (DEBUGLEVEL>5) err_printf("delta[%ld] = %Ps\n",i,p1);
     732             :     /* g = prod (X - delta[i])
     733             :      * if g o h = 0 (pol), we'll have h(Ai[j]) = delta[i] for all j */
     734             :     /* fk[k] belongs to block number whichdelta[k] */
     735      160272 :     for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
     736      160272 :     if (typ(p1) == t_POL) p1 = constant_coeff(p1);
     737      160272 :     d_1_term = addii(d_1_term, p1);
     738             :   }
     739       13930 :   d_1_term = centermod(d_1_term, pe); /* Tr(g) */
     740       13930 :   if (abscmpii(d_1_term, gel(M,3)) > 0) {
     741        8820 :     if (DEBUGLEVEL>1) err_printf("d-1 test failed\n");
     742        8820 :     return NULL;
     743             :   }
     744        5110 :   g = FqV_roots_to_pol(delta, T, pe, 0);
     745        5110 :   g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
     746        5110 :   if (!ok_coeffs(g,M)) {
     747        4494 :     if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
     748        4494 :     if (DEBUGLEVEL>1) err_printf("coeff too big for pol g(x)\n");
     749        4494 :     return NULL;
     750             :   }
     751         616 :   if (!FpX_is_squarefree(g, p)) {
     752          84 :     if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
     753          84 :     if (DEBUGLEVEL>1) err_printf("changing f(x): p divides disc(g)\n");
     754          84 :     compute_data(B);
     755          84 :     return subfield(A, B);
     756             :   }
     757             : 
     758         532 :   lf = lg(firstroot); listdelta = cgetg(lf, t_VEC);
     759         532 :   for (i=1; i<lf; i++) listdelta[i] = delta[whichdelta[firstroot[i]]];
     760         532 :   if (DEBUGLEVEL) err_printf("candidate = %Ps\n", g);
     761         532 :   e = embedding(g, B->DATA, B->S, B->PD->den, listdelta);
     762         532 :   if (!e) return NULL;
     763         532 :   if (DEBUGLEVEL) err_printf("... OK!\n");
     764         532 :   return _subfield(g, e);
     765             : }
     766             : 
     767             : /* L list of current subfields, test whether potential block D is a block,
     768             :  * if so, append corresponding subfield */
     769             : static GEN
     770       13846 : test_block(blockdata *B, GEN L, GEN D)
     771             : {
     772       13846 :   pari_sp av = avma;
     773       13846 :   GEN sub = subfield(D, B);
     774       13846 :   if (sub) {
     775         532 :     GEN old = L;
     776         532 :     L = gclone( L? shallowconcat(L, sub): sub );
     777         532 :     guncloneNULL(old);
     778             :   }
     779       13846 :   return gc_const(av,L);
     780             : }
     781             : 
     782             : /* subfields of degree d */
     783             : static GEN
     784         455 : subfields_of_given_degree(blockdata *B)
     785             : {
     786         455 :   pari_sp av = avma;
     787             :   GEN L;
     788             : 
     789         455 :   if (DEBUGLEVEL) err_printf("\n* Look for subfields of degree %ld\n\n", B->d);
     790         455 :   B->DATA = NULL; compute_data(B);
     791         455 :   L = calc_block(B, B->S->Z, cgetg(1,t_VEC), NULL);
     792         455 :   if (DEBUGLEVEL>9)
     793           0 :     err_printf("\nSubfields of degree %ld: %Ps\n", B->d, L? L: cgetg(1,t_VEC));
     794         455 :   clone_unlock(B->DATA); return gc_const(av,L);
     795             : }
     796             : 
     797             : static GEN
     798         266 : fix_var(GEN x, long v)
     799             : {
     800         266 :   long i, l = lg(x);
     801         266 :   if (!v) return x;
     802           0 :   for (i=1; i<l; i++) { GEN t = gel(x,i); setvarn(gel(t,1),v); setvarn(gel(t,2),v); }
     803           0 :   return x;
     804             : }
     805             : 
     806             : static void
     807         266 : subfields_poldata(GEN nf, GEN T, poldata *PD)
     808             : {
     809             :   GEN L, dis;
     810             : 
     811         266 :   PD->pol = T;
     812         266 :   if (nf)
     813             :   {
     814         161 :     PD->den = nf_get_zkden(nf);
     815         161 :     PD->roo = nf_get_roots(nf);
     816         161 :     PD->dis = mulii(absi_shallow(nf_get_disc(nf)), sqri(nf_get_index(nf)));
     817             :   }
     818             :   else
     819             :   {
     820         105 :     PD->den = initgaloisborne(T,NULL,nbits2prec(bit_accuracy(ZX_max_lg(T))), &L,NULL,&dis);
     821         105 :     PD->roo = L;
     822         105 :     PD->dis = absi_shallow(dis);
     823             :   }
     824         266 : }
     825             : 
     826             : static GEN nfsubfields_fa(GEN nf, long d);
     827             : static GEN
     828         483 : subfieldsall(GEN nf0)
     829             : {
     830         483 :   pari_sp av = avma;
     831             :   long N, ld, i, v;
     832             :   GEN nf, G, T, dg, LSB, NLSB;
     833             :   poldata PD;
     834             :   primedata S;
     835             :   blockdata B;
     836             : 
     837             :   /* much easier if nf is Galois (WSS) */
     838         483 :   G = galoisinit(nf0, NULL);
     839         483 :   T = get_nfpol(nf0, &nf);
     840         483 :   if (G != gen_0)
     841             :   {
     842             :     GEN L, S;
     843             :     long l;
     844          28 :     L = lift_shallow( galoissubfields(G, 0, varn(T)) );
     845          28 :     l = lg(L); S = cgetg(l, t_VECSMALL);
     846          28 :     for (i=1; i<l; i++) S[i] = lg(gmael(L,i,1));
     847          28 :     return gerepilecopy(av, vecpermute(L, vecsmall_indexsort(S)));
     848             :   }
     849         455 :   v = varn(T); N = degpol(T);
     850         455 :   dg = divisorsu(N); ld = lg(dg)-1;
     851         455 :   LSB = _subfield(pol_x(v), pol_0(v));
     852         455 :   if (ld <= 2)
     853             :   {
     854         189 :     if (ld == 2) LSB = shallowconcat(LSB, _subfield(T, pol_x(v)));
     855         189 :     return gerepilecopy(av, LSB);
     856             :   }
     857         266 :   if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
     858         266 :   if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, 0); }
     859         259 :   subfields_poldata(nf, T, &PD);
     860             : 
     861         259 :   if (DEBUGLEVEL) err_printf("\n***** Entering subfields\n\npol = %Ps\n",T);
     862         259 :   B.PD = &PD;
     863         259 :   B.S  = &S;
     864         259 :   B.N  = N;
     865         707 :   for (i=ld-1; i>1; i--)
     866             :   {
     867         448 :     B.size  = dg[i];
     868         448 :     B.d = N / B.size;
     869         448 :     NLSB = subfields_of_given_degree(&B);
     870         448 :     if (NLSB) { LSB = gconcat(LSB, NLSB); gunclone(NLSB); }
     871             :   }
     872         259 :   (void)delete_var(); /* from init_primedata() */
     873         259 :   LSB = shallowconcat(LSB, _subfield(T, pol_x(0)));
     874         259 :   if (DEBUGLEVEL) err_printf("\n***** Leaving subfields\n\n");
     875         259 :   return fix_var(gerepilecopy(av, LSB), v);
     876             : }
     877             : 
     878             : GEN
     879         742 : nfsubfields(GEN nf0, long d)
     880             : {
     881         742 :   pari_sp av = avma;
     882             :   long N, v0;
     883             :   GEN nf, LSB, T, G;
     884             :   poldata PD;
     885             :   primedata S;
     886             :   blockdata B;
     887             : 
     888         742 :   if (typ(nf0)==t_VEC && lg(nf0)==3) return nfsubfields_fa(nf0, d);
     889         532 :   if (!d) return subfieldsall(nf0);
     890             : 
     891             :   /* treat trivial cases */
     892          49 :   T = get_nfpol(nf0, &nf); v0 = varn(T); N = degpol(T);
     893          49 :   RgX_check_ZX(T,"nfsubfields");
     894          49 :   if (d == N) return gerepilecopy(av, _subfield(T, pol_x(v0)));
     895          35 :   if (d == 1) return gerepilecopy(av, _subfield(pol_x(v0), zeropol(v0)));
     896          21 :   if (d < 1 || d > N || N % d) return cgetg(1,t_VEC);
     897             : 
     898             :   /* much easier if nf is Galois (WSS) */
     899          14 :   G = galoisinit(nf0, NULL);
     900          14 :   if (G != gen_0)
     901             :   { /* Bingo */
     902           7 :     GEN L = galoissubgroups(G), F;
     903           7 :     long k,i, l = lg(L), o = N/d;
     904           7 :     F = cgetg(l, t_VEC);
     905           7 :     k = 1;
     906          49 :     for (i=1; i<l; i++)
     907             :     {
     908          42 :       GEN H = gel(L,i);
     909          42 :       if (group_order(H) == o)
     910          21 :         gel(F,k++) = lift_shallow(galoisfixedfield(G, gel(H,1), 0, v0));
     911             :     }
     912           7 :     setlg(F, k);
     913           7 :     return gerepilecopy(av, F);
     914             :   }
     915           7 :   if (varn(T) != 0) { T = leafcopy(T); setvarn(T, 0); }
     916           7 :   if (!choose_prime(&S, T)) { set_avma(av); return nfsubfields_fa(nf0, d); }
     917           7 :   subfields_poldata(nf, T, &PD);
     918           7 :   B.PD = &PD;
     919           7 :   B.S  = &S;
     920           7 :   B.N  = N;
     921           7 :   B.d  = d;
     922           7 :   B.size = N/d;
     923           7 :   LSB = subfields_of_given_degree(&B);
     924           7 :   (void)delete_var(); /* from init_primedata */
     925           7 :   set_avma(av);
     926           7 :   if (!LSB) return cgetg(1, t_VEC);
     927           7 :   G = gcopy(LSB); gunclone(LSB);
     928           7 :   return fix_var(G, v0);
     929             : }
     930             : 
     931             : /******************************/
     932             : /*                            */
     933             : /*    Maximal CM subfield     */
     934             : /*     Aurel Page (2019)      */
     935             : /*                            */
     936             : /******************************/
     937             : 
     938             : /* ero: maximum exponent+1 of roots of pol */
     939             : static GEN
     940        2695 : try_subfield_generator(GEN pol, GEN v, long e, long p, long ero)
     941             : {
     942             :   GEN a, P, Q;
     943             :   long d, bound, i, B, bi, ed;
     944             : 
     945        2695 :   a = gtopolyrev(v, varn(pol));
     946        2695 :   P = Flxq_charpoly(ZX_to_Flx(a,p), ZX_to_Flx(pol,p), p);
     947        2695 :   Flx_ispower(P, e, p, &Q);
     948        2695 :   if (!Flx_is_squarefree(Q,p)) return NULL;
     949        1337 :   d = degree(pol)/e;
     950        1337 :   B = 0;
     951       21938 :   for (i=1; i<lg(v); i++)
     952             :   {
     953       20601 :     bi = (i-1)*ero + expi(gel(v,i));
     954       20601 :     if (bi > B) B = bi;
     955             :   }
     956        1337 :   ed = expu(d);
     957        1337 :   B += ed+1;
     958        1337 :   bound = 0;
     959        6741 :   for (i=0; 2*i<=d; i++)
     960             :   {
     961        5404 :     if (!i) bi = d*B;
     962        4067 :     else    bi = (d-i)*B + i*(3+ed-expu(i));
     963        5404 :     if (bi > bound) bound = bi;
     964             :   }
     965        1337 :   Q = ZXQ_minpoly(a,pol,d,bound);
     966        1337 :   return mkvec2(Q, a);
     967             : }
     968             : 
     969             : /* subfield sub of nf of degree d assuming:
     970             :    - V is contained in sub
     971             :    - V is not contained in a proper subfield of sub
     972             :    ero: maximum exponent+1 of roots of pol
     973             :    output as nfsubfields:
     974             :    - pair [g,h] where g absolute equation for the  subfield and h expresses
     975             :    - one of the roots of g in terms of the generator of nf
     976             : */
     977             : static GEN
     978        1491 : subfield_generator(GEN pol, GEN V, long d, long ero)
     979             : {
     980             :   long p, i, e;
     981        1491 :   GEN a = NULL, v = cgetg(lg(V),t_COL), B;
     982             : 
     983        1491 :   if (d==1) return mkvec2(pol_x(varn(pol)), gen_0);
     984        1337 :   e = degree(pol)/d;
     985        1337 :   p = 1009;
     986        2688 :   for (i=1; i<lg(V); i++)
     987             :   {
     988        2667 :     a = try_subfield_generator(pol, gel(V,i), e, p, ero);
     989        2667 :     if (a) return a;
     990        1351 :     p = unextprime(p+1);
     991             :   }
     992          21 :   B = stoi(10);
     993             :   while(1)
     994             :   {
     995          35 :     for (i=1; i<lg(v); i++) gel(v,i) = randomi(B);
     996          28 :     a = try_subfield_generator(pol, QM_QC_mul(V,v), e, p, ero);
     997          28 :     if (a) return a;
     998           7 :     p = unextprime(p+1);
     999             :   }
    1000             :   return NULL; /* not reached */
    1001             : }
    1002             : 
    1003             : static GEN
    1004       29442 : polmod_to_col(GEN P, long dx, long dy)
    1005             : {
    1006             :   GEN res, c;
    1007             :   long i, j, k;
    1008             : 
    1009       29442 :   P = liftpol(P);
    1010       29442 :   res = zerocol((dx+1)*(dy+1));
    1011       29442 :   k = 1;
    1012      110474 :   for(i=0; i<=dy; i++)
    1013             :   {
    1014       81032 :     c = polcoef_i(P,i,-1);
    1015     1907444 :     for (j=0; j<=dx; j++)
    1016             :     {
    1017     1826412 :       gel(res,k) = polcoef_i(c,j,-1);
    1018     1826412 :       k++;
    1019             :     }
    1020             :   }
    1021       29442 :   return res;
    1022             : }
    1023             : 
    1024             : /* lambda: t_VEC of t_INT; 0 means ignore this factor */
    1025             : static GEN
    1026        1883 : twoembequation(GEN pol, GEN fa, GEN lambda)
    1027             : {
    1028             :   GEN m, vpolx, vpoly, x, y, C, modpol;
    1029             :   long i,j, dx, dy;
    1030             : 
    1031        1883 :   x = pol_x(varn(pol));
    1032        1883 :   dx = degree(pol);
    1033        1883 :   y = pol_x(varn(gel(fa,1)));
    1034        1883 :   modpol = mkpolmod(gen_1,pol);
    1035             : 
    1036        1883 :   lambda = shallowcopy(lambda);
    1037        1883 :   fa = shallowcopy(fa);
    1038        1883 :   j = 1;
    1039       21623 :   for (i=1; i<lg(fa); i++)
    1040       19740 :     if (signe(gel(lambda,i)))
    1041             :     {
    1042        1904 :       gel(lambda,j) = negi(gel(lambda,i));
    1043        1904 :       gel(fa,j) = gmul(gel(fa,i),modpol);
    1044        1904 :       j++;
    1045             :     }
    1046        1883 :   setlg(lambda, j);
    1047        1883 :   setlg(fa, j);
    1048             : 
    1049        1883 :   vpolx = cgetg(lg(fa),t_VEC);
    1050        3787 :   for (i=1; i<lg(fa); i++)
    1051        1904 :     gel(vpolx,i) = gmul(modpol,mkpolmod(gen_1,gel(fa,i)));
    1052        1883 :   vpoly = gcopy(vpolx);
    1053             : 
    1054        1883 :   m = cgetg(degree(pol)+1,t_MAT);
    1055       31171 :   for (j=1; j<lg(m); j++)
    1056             :   {
    1057       29288 :     C = zerovec(lg(fa)-1);
    1058       58730 :     for(i=1; i<lg(fa); i++)
    1059             :     {
    1060       29442 :       dy = degree(gel(fa,i));
    1061       29442 :       gel(C,i) = polmod_to_col(gadd(gel(vpolx,i),gmul(gel(lambda,i),gel(vpoly,i))), dx, dy);
    1062       29442 :       gel(vpolx,i) = gmul(gel(vpolx,i),x);
    1063       29442 :       gel(vpoly,i) = gmul(gel(vpoly,i),y);
    1064             :     }
    1065       29288 :     gel(m,j) = shallowconcat1(C);
    1066             :   }
    1067        1883 :   return QM_ker(m);
    1068             : }
    1069             : 
    1070             : static void
    1071         588 : subfields_cleanup(GEN* nf, GEN* pol, long* n, GEN* fa)
    1072             : {
    1073         588 :   *fa = NULL;
    1074         588 :   if (typ(*nf) != t_VEC && typ(*nf) != t_POL) pari_err_TYPE("subfields_cleanup", *nf);
    1075         581 :   if (typ(*nf) == t_VEC && lg(*nf) == 3)
    1076             :   {
    1077         266 :     *fa = gel(*nf,2);
    1078         266 :     *nf = gel(*nf,1);
    1079         266 :     if (typ(*fa)!=t_MAT || lg(*fa)!=3)
    1080          14 :       pari_err_TYPE("subfields_cleanup [fa should be a factorisation matrix]", *fa);
    1081             :   }
    1082         567 :   if (typ(*nf) == t_POL)
    1083             :   {
    1084         434 :     *pol = *nf;
    1085         434 :     *nf = NULL;
    1086         434 :     if (!gequal1(leading_coeff(*pol))) pari_err_TYPE("subfields_cleanup [not monic]", *pol);
    1087         427 :     *n = degree(*pol);
    1088         427 :     if (*n<=0) pari_err_TYPE("subfields_cleanup [constant polynomial]", *pol);
    1089             :   }
    1090             :   else
    1091             :   {
    1092         133 :     *nf = checknf(*nf);
    1093         112 :     *pol = nf_get_pol(*nf);
    1094         112 :     *n = degree(*pol);
    1095             :   }
    1096         532 : }
    1097             : 
    1098             : static GEN
    1099         134 : rootsuptoconj(GEN pol, long prec)
    1100             : {
    1101             :   GEN ro;
    1102             :   long n, i;
    1103         134 :   ro = roots(pol,prec);
    1104         134 :   n = lg(ro)-1;
    1105         917 :   for (i=1; i<=n/2; i++)
    1106         783 :     gel(ro,i) = gel(ro,2*i-1);
    1107         134 :   setlg(ro,n/2+1);
    1108         134 :   return ro;
    1109             : }
    1110             : static GEN
    1111         238 : cmsubfield_get_roots(GEN pol, GEN nf, long n, long* r2, long *prec)
    1112             : {
    1113             :   GEN ro;
    1114         238 :   if (nf)
    1115             :   {
    1116         112 :     if (nf_get_r1(nf)) return NULL;
    1117         112 :     *r2 = nf_get_r2(nf);
    1118         112 :     *prec = nf_get_prec(nf);
    1119         112 :     ro = nf_get_roots(nf);
    1120             :   }
    1121             :   else
    1122             :   {
    1123         126 :     if (n%2 || sturm(pol)) return NULL;
    1124          91 :     *r2 = n/2;
    1125          91 :     *prec = MEDDEFAULTPREC;
    1126          91 :     ro = rootsuptoconj(pol, *prec);
    1127             :   }
    1128         203 :   return ro;
    1129             : }
    1130             : 
    1131             : static GEN
    1132         357 : subfields_get_fa(GEN pol, GEN nf, GEN fa)
    1133             : {
    1134             :   GEN poly;
    1135         357 :   long var = fetch_var_higher();
    1136         357 :   if (!fa)
    1137             :   {
    1138         168 :     poly = shallowcopy(pol);
    1139         168 :     setvarn(poly,var);
    1140         168 :     if (nf) fa = nffactor(nf, poly);
    1141         105 :     else    fa = nffactor(pol, poly);
    1142             :   }
    1143         357 :   return liftpol(gel(fa,1));
    1144             : }
    1145             : 
    1146             : static long
    1147         217 : subfields_get_ero(GEN pol, GEN nf)
    1148             : {
    1149             :   GEN ro;
    1150             :   long ero, e, i;
    1151         217 :   if (nf)   ro = nf_get_roots(nf);
    1152         217 :   else      ro = roots(pol, LOWDEFAULTPREC);
    1153         217 :   ero = 0;
    1154        2933 :   for (i=1; i<lg(ro); i++)
    1155             :   {
    1156        2716 :     e = 1+gexpo(gel(ro,i));
    1157        2716 :     if (e > ero) ero = e;
    1158             :   }
    1159         217 :   return ero+1;
    1160             : };
    1161             : 
    1162             : static GEN
    1163          49 : try_imag(GEN x, GEN c, GEN pol, long v, ulong p, GEN emb, GEN galpol)
    1164             : {
    1165             :   GEN a;
    1166          49 :   a = Q_primpart(RgX_sub(RgX_RgXQ_eval(x,c,pol),x));
    1167          49 :   if (Flx_is_squarefree(Flxq_charpoly(ZX_to_Flx(a,p),ZX_to_Flx(pol,p),p),p))
    1168             :   {
    1169          21 :     emb = RgX_RgXQ_eval(a, emb, galpol);
    1170          21 :     pol = ZXQ_charpoly(a, pol, v);
    1171          21 :     return mkvec2(pol, emb);
    1172             :   }
    1173          28 :   return NULL;
    1174             : }
    1175             : 
    1176             : static GEN
    1177          35 : galoissubfieldcm(GEN G)
    1178             : {
    1179          35 :   pari_sp av = avma;
    1180             :   GEN c, H, elts, g, Hset, c2, gene, sub, pol, M, emb, a, galpol, B, b;
    1181             :   long n, i, j, nH, ind, v, d, vsub;
    1182          35 :   ulong p = 1009;
    1183             : 
    1184          35 :   galpol = gal_get_pol(G);
    1185          35 :   n = degree(galpol);
    1186          35 :   v = varn(galpol);
    1187          35 :   c = galois_get_conj(G);
    1188             :   /* compute the list of c*g*c*g^(-1) : product of all pairs of conjugations
    1189             :    * maximal CM subfield is the field fixed by those elements, if c does not
    1190             :    * belong to the group they generate */
    1191          35 :   checkgroup(G, &elts);
    1192          35 :   elts = gen_sort(elts,(void*)vecsmall_lexcmp,cmp_nodata);
    1193          35 :   H = vecsmall_ei(n,1); /* indices of elements of H */
    1194          35 :   Hset = zero_F2v(n);
    1195          35 :   F2v_set(Hset,1);
    1196          35 :   nH = 1;
    1197         168 :   for (i=2; i<=n; i++)
    1198             :   {
    1199         133 :     g = gel(elts,i);
    1200         133 :     c2 = perm_mul(c,perm_conj(g,c));
    1201         133 :     if (!F2v_coeff(Hset,c2[1]))
    1202             :     {
    1203          14 :       nH++;
    1204          14 :       H[nH] = c2[1];
    1205          14 :       F2v_set(Hset,c2[1]);
    1206             :     }
    1207             :   }
    1208             :   /* group generated */
    1209          35 :   gene = gcopy(H);
    1210          35 :   setlg(gene,nH+1);
    1211          35 :   i = 1; /* last element that has been multiplied by the generators */
    1212          84 :   while (i < nH)
    1213             :   {
    1214          56 :     for (j=1; j<lg(gene); j++)
    1215             :     {
    1216          42 :       g = gel(elts,gene[j]);
    1217          42 :       ind = g[H[i]]; /* index of the product */
    1218          42 :       if (!F2v_coeff(Hset,ind))
    1219             :       {
    1220           0 :         nH++;
    1221           0 :         if (ind==c[1] || 2*nH>n) { set_avma(av); return gen_0; }
    1222           0 :         H[nH] = ind;
    1223           0 :         F2v_set(Hset,ind);
    1224             :       }
    1225             :     }
    1226          14 :     i++;
    1227             :   }
    1228          35 :   H = cgetg(lg(gene), t_VEC);
    1229          84 :   for (i=1; i<lg(H); i++)
    1230          49 :     gel(H,i) = gel(elts,gene[i]);
    1231          35 :   vsub = fetch_var();
    1232          35 :   sub = galoisfixedfield(G, H, 2, vsub);
    1233             : 
    1234             :   /* compute a totally imaginary generator */
    1235          35 :   pol = gel(sub,1);
    1236          35 :   emb = liftpol(gel(sub,2));
    1237          35 :   d = degree(pol);
    1238          35 :   if (!(ZX_deflate_order(pol)%2) && sturm(RgX_deflate(pol,2))==d/2)
    1239             :   {
    1240          14 :     setvarn(pol,v);
    1241          14 :     delete_var();
    1242          14 :     return mkvec2(pol,emb);
    1243             :   }
    1244             : 
    1245             :   /* compute action of c on the subfield from that on the large field */
    1246          21 :   c = galoispermtopol(G,c);
    1247          21 :   if (d<n)
    1248             :   {
    1249           0 :     c = RgX_to_RgC(c,n);
    1250           0 :     M = cgetg(d+1,t_MAT);
    1251           0 :     gel(M,1) = col_ei(n,1);
    1252           0 :     a = pol_1(v);
    1253           0 :     for (i=2; i<=d; i++)
    1254             :     {
    1255           0 :       a = ZX_rem(ZX_mul(a,emb), galpol);
    1256           0 :       gel(M,i) = RgX_to_RgC(a,n);
    1257             :     }
    1258           0 :     c = RgM_RgC_invimage(M, c);
    1259           0 :     c = RgV_to_RgX(c, vsub);
    1260             :   }
    1261          21 :   else setvarn(c,vsub);
    1262             : 
    1263             :   /* search for a generator of the form c(b)-b */
    1264          42 :   for (i=1; i<d; i++)
    1265             :   {
    1266          35 :     a = try_imag(pol_xn(i,vsub),c,pol,v,p,emb,galpol);
    1267          35 :     if (a) { delete_var(); return a; }
    1268          21 :     p = unextprime(p+1);
    1269             :   }
    1270           7 :   B = stoi(10);
    1271           7 :   b = pol_xn(d-1,vsub);
    1272             :   while(1)
    1273             :   {
    1274          21 :     for (i=2; i<lg(b); i++) gel(b,i) = randomi(B);
    1275          14 :     a = try_imag(b,c,pol,v,p,emb,galpol);
    1276          14 :     if (a) { delete_var(); return a; }
    1277           7 :     p = unextprime(p+1);
    1278             :   }
    1279             :   return NULL; /* not reached */
    1280             : }
    1281             : 
    1282             : static GEN
    1283          28 : quadsubfieldcm(GEN pol)
    1284             : {
    1285          28 :   GEN a = polcoef_i(pol,1,-1), b = polcoef_i(pol,0,-1), res;
    1286             :   long i;
    1287          28 :   if (smodis(a,2))
    1288          14 :     res = mkvec2(
    1289             :           mkpoln(3,gen_1,gen_0,subii(mulsi(4,b),sqri(a))),
    1290             :           mkpoln(2,gen_2,a)
    1291             :           );
    1292             :   else
    1293             :   {
    1294          14 :     a = divis(a,2);
    1295          14 :     res = mkvec2(
    1296             :           mkpoln(3,gen_1,gen_0,subii(b,sqri(a))),
    1297             :           mkpoln(2,gen_1,a)
    1298             :           );
    1299             :   }
    1300          28 :   for (i=1; i<=2; i++) setvarn(gel(res,i),varn(pol));
    1301          28 :   return res;
    1302             : }
    1303             : 
    1304             : GEN
    1305         294 : nfsubfieldscm(GEN nf)
    1306             : {
    1307         294 :   pari_sp av = avma;
    1308             :   GEN fa, lambda, V, res, ro, a, aa, ev, minev, pol, G;
    1309         294 :   long i, j, n, r2, minj=0, prec, emax, emin, e, precbound, ero;
    1310             : 
    1311         294 :   subfields_cleanup(&nf, &pol, &n, &fa);
    1312         238 :   ro = cmsubfield_get_roots(pol, nf, n, &r2, &prec);
    1313         238 :   if (!ro) { set_avma(av); return gen_0; }
    1314             :   /* now r2 == 2*n */
    1315             : 
    1316         203 :   if (n==2) return gerepilecopy(av, quadsubfieldcm(pol));
    1317         175 :   G = galoisinit(nf? nf: pol, NULL);
    1318         175 :   if (G != gen_0) return gerepilecopy(av, galoissubfieldcm(G));
    1319             : 
    1320         140 :   ero = 0;
    1321         819 :   for (i=1; i<lg(ro); i++)
    1322             :   {
    1323         679 :     e = 1+gexpo(gel(ro,i));
    1324         679 :     if (e > ero) ero = e;
    1325             :   }
    1326         140 :   ero++;
    1327         140 :   fa = subfields_get_fa(pol, nf, fa);
    1328             : 
    1329         140 :   emax = 1;
    1330         140 :   emin = -1;
    1331         819 :   for (i=1; i<lg(ro); i++)
    1332        2723 :     for (j=i+1; j<lg(ro); j++)
    1333             :     {
    1334        2044 :       e = gexpo(gsub(gel(ro,i),gel(ro,j)));
    1335        2044 :       if (e > emax) emax = e;
    1336        2044 :       if (e < emin) emin = e;
    1337             :     }
    1338         140 :   precbound = n*(emax-emin) + gexpo(fa) + n*n + 5;
    1339         140 :   precbound = 3 + precbound/BITS_IN_LONG;
    1340         140 :   if (prec < precbound)
    1341             :   {
    1342          43 :     prec = precbound;
    1343          43 :     ro = rootsuptoconj(pol, prec);
    1344             :   }
    1345             : 
    1346         140 :   lambda = zerovec(lg(fa)-1);
    1347         819 :   for (i=1; i<=r2; i++)
    1348             :   {
    1349         679 :     a = gel(ro,i);
    1350         679 :     aa = conj_i(a);
    1351        5271 :     for (j=1; j<lg(fa); j++)
    1352             :     {
    1353        4592 :       ev = cxnorm(poleval(poleval(gel(fa,j),aa),a));
    1354        4592 :       if (j==1 || cmprr(minev,ev)>0) { minj = j; minev = ev; }
    1355             :     }
    1356         679 :     gel(lambda,minj) = gen_m1;
    1357             :   }
    1358             : 
    1359         140 :   V = twoembequation(pol, fa, lambda);
    1360         140 :   if (lg(V)==1) { set_avma(av); delete_var(); return gen_0; }
    1361         133 :   res = subfield_generator(pol, V, 2*(lg(V)-1), ero);
    1362         133 :   delete_var();
    1363         133 :   return gerepilecopy(av, res);
    1364             : }
    1365             : 
    1366             : static int
    1367       17962 : field_is_contained(GEN V, GEN W, int strict)
    1368             : {
    1369             :   GEN VW;
    1370       17962 :   ulong p = 1073741827;
    1371             :   /* distinct overfield must have different dimension */
    1372       17962 :   if (strict && lg(V) == lg(W)) return 0;
    1373             :   /* dimension of overfield must be multiple */
    1374       13538 :   if ((lg(W)-1) % (lg(V)-1)) return 0;
    1375        9709 :   VW = shallowconcat(V,W);
    1376        9709 :   if (Flm_rank(ZM_to_Flm(VW,p),p) > lg(W)-1) return 0;
    1377        3836 :   return ZM_rank(VW) == lg(W)-1;
    1378             : }
    1379             : 
    1380             : /***********************************************/
    1381             : /*                                             */
    1382             : /*    Maximal, generating, all subfields       */
    1383             : /*             Aurel Page (2019)               */
    1384             : /*     after van Hoeij, Klueners, Novocin      */
    1385             : /*  Journal of Symbolic Computation 52 (2013)  */
    1386             : /*                                             */
    1387             : /***********************************************/
    1388             : 
    1389             : const long subf_MAXIMAL = 1; /* return the maximal subfields */
    1390             : const long subf_GENERATING = 2; /* return the generating subfields */
    1391             : static GEN
    1392         217 : maxgen_subfields(GEN pol, GEN fa, long flag)
    1393             : {
    1394         217 :   pari_sp av = avma;
    1395         217 :   GEN principal, ismax, isgene, Lmax = NULL, Lgene, res, V, W, W1;
    1396         217 :   long i, i2, j, flmax, flgene, nbmax = 0, nbgene = 0, dup;
    1397             : 
    1398         217 :   if (!flag) return cgetg(1,t_VEC);
    1399         217 :   flmax = (flag & subf_MAXIMAL)!=0;
    1400         217 :   flgene = (flag & subf_GENERATING)!=0;
    1401             : 
    1402             :   /* compute principal subfields */
    1403         217 :   principal = cgetg(lg(fa),t_VEC);
    1404        1960 :   for (i=1; i<lg(fa); i++)
    1405        1743 :     gel(principal,i) = twoembequation(pol, fa, vec_ei(lg(fa)-1,i));
    1406         217 :   principal = gen_sort_uniq(principal, (void*)&cmp_universal, &cmp_nodata);
    1407             :   /* remove nf and duplicates (sort_uniq possibly not enough) */
    1408         217 :   i2 = 1;
    1409        1309 :   for (i=1; i<lg(principal)-1; i++)
    1410             :   {
    1411        1092 :     V = gel(principal,i);
    1412        1092 :     dup = 0;
    1413        1092 :     j = i2-1;
    1414        3493 :     while(j>0 && lg(gel(principal,j))==lg(V))
    1415             :     {
    1416        1309 :       if (field_is_contained(gel(principal,j),V,0)) { dup=1; break; }
    1417        1309 :       j--;
    1418             :     }
    1419        1092 :     if (!dup)
    1420             :     {
    1421        1092 :       gel(principal,i2) = V;
    1422        1092 :       i2++;
    1423             :     }
    1424             :   }
    1425         217 :   setlg(principal, i2);
    1426             : 
    1427             :   /* a subfield is generating iff all overfields contain the first overfield */
    1428         217 :   ismax = cgetg(lg(principal),t_VECSMALL);
    1429         217 :   isgene = cgetg(lg(principal),t_VECSMALL);
    1430        1309 :   for (i=1; i<lg(principal); i++)
    1431             :   {
    1432        1092 :     V = gel(principal,i);
    1433        1092 :     ismax[i] = flmax;
    1434        1092 :     isgene[i] = flgene;
    1435        1092 :     W1 = NULL; /* intersection of strict overfields */
    1436        4137 :     for (j=i+1; j<lg(principal); j++)
    1437             :     {
    1438        3185 :       W = gel(principal,j);
    1439        3185 :       if (!field_is_contained(V,W,1)) continue;
    1440         567 :       ismax[i] = 0;
    1441         567 :       if (!flgene) break;
    1442         462 :       if (!W1) { W1 = W; continue; }
    1443         189 :       if (!field_is_contained(W1,W,1))
    1444             :       {
    1445          63 :         W1 = intersect(W1,W);
    1446          63 :         if (lg(W1)==lg(V)) { isgene[i]=0; break; }
    1447             :       }
    1448             :     }
    1449             :   }
    1450             : 
    1451        1309 :   for (i=1; i<lg(principal); i++)
    1452             :   {
    1453        1092 :     nbmax += ismax[i];
    1454        1092 :     nbgene += isgene[i];
    1455             :   }
    1456             : 
    1457         217 :   if (flmax)
    1458             :   {
    1459          49 :     Lmax = cgetg(nbmax+1, t_VEC);
    1460          49 :     j=1;
    1461         259 :     for (i=1; i<lg(principal); i++)
    1462         210 :       if (ismax[i])
    1463             :       {
    1464         105 :         gel(Lmax,j) = gel(principal,i);
    1465         105 :         j++;
    1466             :       }
    1467             :   }
    1468             : 
    1469         217 :   if (flgene)
    1470             :   {
    1471         168 :     Lgene = cgetg(nbgene+1, t_VEC);
    1472         168 :     j=1;
    1473        1050 :     for (i=1; i<lg(principal); i++)
    1474         882 :       if (isgene[i])
    1475             :       {
    1476         847 :         gel(Lgene,j) = gel(principal,i);
    1477         847 :         j++;
    1478             :       }
    1479             :   }
    1480             : 
    1481         217 :   if (!flgene) res = Lmax;
    1482         168 :   else if (!flmax) res = Lgene;
    1483           0 :   else res = mkvec2(Lmax,Lgene);
    1484         217 :   return gerepilecopy(av, res);
    1485             : }
    1486             : 
    1487             : GEN
    1488          77 : nfsubfieldsmax(GEN nf)
    1489             : {
    1490          77 :   pari_sp av = avma;
    1491             :   GEN pol, fa, Lmax, V;
    1492             :   long n, i, ero;
    1493             : 
    1494          77 :   subfields_cleanup(&nf, &pol, &n, &fa);
    1495          77 :   if (n==1) { set_avma(av); return cgetg(1,t_VEC); }
    1496          70 :   if (uisprime(n)) return gerepilecopy(av,mkvec(mkvec2(pol_x(varn(pol)),gen_0)));
    1497          49 :   ero = subfields_get_ero(pol, nf);
    1498          49 :   fa = subfields_get_fa(pol, nf, fa);
    1499          49 :   Lmax = maxgen_subfields(pol, fa, subf_MAXIMAL);
    1500         154 :   for (i=1; i<lg(Lmax); i++)
    1501             :   {
    1502         105 :     V = gel(Lmax,i);
    1503         105 :     gel(Lmax,i) = subfield_generator(pol, V, lg(V)-1, ero);
    1504             :   }
    1505          49 :   delete_var();
    1506          49 :   return gerepilecopy(av, Lmax);
    1507             : }
    1508             : 
    1509             : static void
    1510        1575 : heap_climb(GEN* H, long i)
    1511             : {
    1512             :   long j;
    1513        1575 :   if (i==1) return;
    1514        1169 :   j = i/2;
    1515        1169 :   if (cmp_universal(gel(*H,i),gel(*H,j)) > 0)
    1516             :   {
    1517         476 :     swap(gel(*H,i), gel(*H,j));
    1518         476 :     return heap_climb(H,j);
    1519             :   }
    1520             : }
    1521             : 
    1522             : static void
    1523        1099 : heap_push(GEN* H, long *len, GEN x)
    1524             : {
    1525        1099 :   if (*len+1 == lg(*H))
    1526             :   {
    1527          14 :     GEN H2 = zerovec(2*(*len));
    1528             :     long i;
    1529         154 :     for(i=1; i<lg(*H); i++)
    1530         140 :       gel(H2,i) = gel(*H,i);
    1531          14 :     *H = H2;
    1532             :   }
    1533        1099 :   (*len)++;
    1534        1099 :   gel(*H,*len) = x;
    1535        1099 :   return heap_climb(H,*len);
    1536             : }
    1537             : 
    1538             : static void
    1539        2324 : heap_descend(GEN* H, long len, long i)
    1540             : {
    1541        2324 :   long maxi = i, j = 2*i;
    1542        2324 :   if (j > len) return;
    1543        1232 :   if (cmp_universal(gel(*H,j),gel(*H,i)) > 0) maxi = j;
    1544        1232 :   j++;
    1545        1232 :   if (j<=len && cmp_universal(gel(*H,j),gel(*H,maxi))>0) maxi = j;
    1546        1232 :   if (maxi == i) return;
    1547        1057 :   swap(gel(*H,i), gel(*H,maxi));
    1548        1057 :   return heap_descend(H,len,maxi);
    1549             : }
    1550             : 
    1551             : static void
    1552        1267 : heap_pop(GEN *H, long *len, GEN* top)
    1553             : {
    1554        1267 :   *top = gel(*H,1);
    1555        1267 :   gel(*H,1) = gel(*H,*len);
    1556        1267 :   (*len)--;
    1557        1267 :   return heap_descend(H,*len,1);
    1558             : };
    1559             : 
    1560             : static GEN
    1561         217 : nfsubfields_fa(GEN nf, long d)
    1562             : {
    1563         217 :   pari_sp av = avma;
    1564             :   GEN pol, fa, gene, res, res2, H, V, v, W, w, data;
    1565             :   long n, r, i, j, nres, len, s, newfield, ero;
    1566             : 
    1567         217 :   subfields_cleanup(&nf, &pol, &n, &fa);
    1568         217 :   if (d && (d<1 || d>n || n%d)) return gerepilecopy(av, cgetg(1,t_VEC));
    1569         238 :   if (!d && uisprime(n)) return gerepilecopy(av, mkvec2(
    1570          14 :         mkvec2(pol_x(varn(pol)),gen_0),
    1571          14 :         mkvec2(pol,pol_x(varn(pol)))
    1572             :         ));
    1573         196 :   if (n==1 || d==1) return gerepilecopy(av,mkvec(mkvec2(pol_x(varn(pol)),gen_0)));
    1574         182 :   if (d==n) return gerepilecopy(av, mkvec(mkvec2(pol,pol_x(varn(pol)))));
    1575         168 :   ero = subfields_get_ero(pol, nf);
    1576         168 :   fa = subfields_get_fa(pol, nf, fa);
    1577         168 :   gene = maxgen_subfields(pol, fa, subf_GENERATING);
    1578             : 
    1579         168 :   if (d)
    1580             :   {
    1581             :     /* keep only generating subfields of degree a multiple of d */
    1582          14 :     j=1;
    1583         147 :     for (i=1; i<lg(gene); i++)
    1584         133 :       if ((lg(gel(gene,i))-1) % d == 0)
    1585             :       {
    1586          98 :         gel(gene,j) = gel(gene,i);
    1587          98 :         j++;
    1588             :       }
    1589          14 :     setlg(gene,j);
    1590             :   }
    1591         168 :   r = lg(gene)-1;
    1592             : 
    1593         168 :   res = zerovec(10);
    1594         168 :   nres = 0;
    1595         168 :   H = zerovec(10);
    1596         168 :   gel(H,1) = mkvec3(matid(n),zero_F2v(r),mkvecsmall(0));
    1597         168 :   len = 1;
    1598             : 
    1599        1603 :   while (len>0)
    1600             :   {
    1601        1267 :     heap_pop(&H, &len, &data);
    1602        1267 :     V = gel(data,1);
    1603        1267 :     v = gel(data,2);
    1604        1267 :     s = gel(data,3)[1];
    1605        5586 :     for (i=s+1; i<=r; i++)
    1606        4319 :       if (!F2v_coeff(v,i))
    1607             :       {
    1608        3360 :         W = vec_Q_primpart(intersect(V, gel(gene,i)));
    1609        3360 :         w = F2v_copy(v);
    1610        3360 :         F2v_set(w, i);
    1611        3360 :         newfield = 1;
    1612       16919 :         for (j=1; j<=r; j++)
    1613       15722 :           if (!F2v_coeff(w,j) && field_is_contained(W,gel(gene,j),1))
    1614             :           {
    1615        3143 :             if (j<i) { newfield = 0; break; }
    1616         980 :             F2v_set(w,j);
    1617             :           }
    1618        3360 :         if (newfield && (!d || (lg(W)-1)%d==0)) heap_push(&H, &len, mkvec3(W,w,mkvecsmall(i)));
    1619             :       }
    1620             : 
    1621        1267 :     if (!d || lg(V)-1==d)
    1622             :     {
    1623        1253 :       nres++;
    1624        1253 :       if (nres == lg(res))
    1625             :       {
    1626          28 :         res2 = zerovec(2*lg(res));
    1627          28 :         for(j=1; j<lg(res); j++) gel(res2,j) = gel(res,j);
    1628          28 :         res = res2;
    1629             :       }
    1630        1253 :       gel(res,nres) = subfield_generator(pol, V, lg(V)-1, ero);
    1631             :     }
    1632             :   }
    1633         168 :   setlg(res,nres+1);
    1634         168 :   vecreverse_inplace(res);
    1635             : 
    1636         168 :   delete_var();
    1637         168 :   return gerepilecopy(av, res);
    1638             : }

Generated by: LCOV version 1.13