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

Generated by: LCOV version 1.14