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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 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.8.0 lcov report (development 16912-212c0f0) Lines: 496 537 92.4 %
Date: 2014-10-20 Functions: 26 26 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 289 356 81.2 %

           Branch data     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 associated 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 associated 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                 :      16555 : calc_block(blockdata *B, GEN Z, GEN Y, GEN SB)
      65                 :            : {
      66                 :      16555 :   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                 :      16555 :   pari_sp av0 = avma;
      69                 :            : 
      70         [ -  + ]:      16555 :   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                 :      16555 :   lnon = minss(BIL, r);
      80                 :      16555 :   e    = new_chunk(BIL);
      81                 :      16555 :   n    = new_chunk(r);
      82                 :      16555 :   non  = new_chunk(lnon);
      83                 :      16555 :   pnon = new_chunk(lnon);
      84                 :      16555 :   pn   = new_chunk(lnon);
      85                 :            : 
      86                 :      16555 :   Zp   = cgetg(lnon,t_VEC);
      87                 :      16555 :   Zpp  = cgetg(lnon,t_VEC); nn = 0;
      88         [ +  + ]:      51870 :   for (i=1; i<r; i++) { n[i] = lg(gel(Z,i))-1; nn += n[i]; }
      89                 :      16555 :   lY = lg(Y); Yp = cgetg(lY+1,t_VEC);
      90         [ +  + ]:      69125 :   for (j=1; j<lY; j++) gel(Yp,j) = gel(Y,j);
      91                 :            : 
      92                 :            :   {
      93                 :      16555 :     pari_sp av = avma;
      94                 :      16555 :     long k = nn / B->size;
      95         [ +  + ]:      36260 :     for (j = 1; j < r; j++)
      96         [ +  + ]:      25620 :       if (n[j] % k) break;
      97         [ +  + ]:      16555 :     if (j == r)
      98                 :            :     {
      99                 :      10640 :       gel(Yp,lY) = Z;
     100                 :      10640 :       SB = print_block_system(B, Yp, SB);
     101                 :      10640 :       avma = av;
     102                 :            :     }
     103                 :            :   }
     104                 :      16555 :   gel(Yp,lY) = Zp;
     105                 :            : 
     106                 :      16555 :   K = divisorsu(n[1]); lK = lg(K);
     107         [ +  + ]:      50113 :   for (i=1; i<lK; i++)
     108                 :            :   {
     109                 :      33558 :     long ngcd = n[1], k = K[i], dk = B->size*k, lpn = 0;
     110         [ +  + ]:      71092 :     for (j=2; j<r; j++)
     111         [ +  + ]:      37534 :       if (n[j]%k == 0)
     112                 :            :       {
     113         [ -  + ]:      37394 :         if (++lpn >= BIL) pari_err_OVERFLOW("calc_block");
     114                 :      37394 :         pn[lpn] = n[j]; pnon[lpn] = j;
     115                 :      37394 :         ngcd = ugcd(ngcd, n[j]);
     116                 :            :       }
     117         [ +  + ]:      33558 :     if (dk % ngcd) continue;
     118                 :      33166 :     T = 1L<<lpn;
     119         [ +  + ]:      33166 :     if (lpn == r-2)
     120                 :            :     {
     121                 :      33096 :       T--; /* done already above --> print_block_system */
     122         [ +  + ]:      33096 :       if (!T) continue;
     123                 :            :     }
     124                 :            : 
     125         [ +  + ]:      20846 :     if (dk == n[1])
     126                 :            :     { /* empty subset, t = 0. Split out for clarity */
     127                 :       8190 :       Zp[1] = Z[1]; setlg(Zp, 2);
     128         [ +  + ]:      21826 :       for (u=1,j=2; j<r; j++) Zpp[u++] = Z[j];
     129                 :       8190 :       setlg(Zpp, u);
     130                 :       8190 :       SB = calc_block(B, Zpp, Yp, SB);
     131                 :            :     }
     132                 :            : 
     133         [ +  + ]:      88193 :     for (t = 1; t < T; t++)
     134                 :            :     { /* loop through all non-empty subsets of [1..lpn] */
     135         [ +  + ]:     354704 :       for (nn=n[1],tp=t, u=1; u<=lpn; u++,tp>>=1)
     136                 :            :       {
     137         [ +  + ]:     287357 :         if (tp&1) { nn += pn[u]; e[u] = 1; } else e[u] = 0;
     138                 :            :       }
     139         [ +  + ]:      67347 :       if (dk != nn) continue;
     140                 :            : 
     141         [ +  + ]:      53725 :       for (j=1; j<r; j++) non[j]=0;
     142                 :       8148 :       Zp[1] = Z[1];
     143         [ +  + ]:      45521 :       for (u=2,j=1; j<=lpn; j++)
     144         [ +  + ]:      37373 :         if (e[j]) { Zp[u] = Z[pnon[j]]; non[pnon[j]] = 1; u++; }
     145                 :       8148 :       setlg(Zp, u);
     146         [ +  + ]:      45577 :       for (u=1,j=2; j<r; j++)
     147         [ +  + ]:      37429 :         if (!non[j]) Zpp[u++] = Z[j];
     148                 :       8148 :       setlg(Zpp, u);
     149                 :       8148 :       SB = calc_block(B, Zpp, Yp, SB);
     150                 :            :     }
     151                 :            :   }
     152                 :      16555 :   avma = av0; return SB;
     153                 :            : }
     154                 :            : 
     155                 :            : /* product of permutations. Put the result in perm1. */
     156                 :            : static void
     157                 :   14475671 : perm_mul_i(GEN perm1, GEN perm2)
     158                 :            : {
     159                 :   14475671 :   long i, N = lg(perm1);
     160                 :   14475671 :   pari_sp av = avma;
     161                 :   14475671 :   GEN perm = new_chunk(N);
     162         [ +  + ]:  405305943 :   for (i=1; i<N; i++) perm[i] = perm1[perm2[i]];
     163         [ +  + ]:  405305943 :   for (i=1; i<N; i++) perm1[i]= perm[i];
     164                 :   14475671 :   avma = av;
     165                 :   14475671 : }
     166                 :            : 
     167                 :            : /* cy is a cycle; compute cy^l as a permutation */
     168                 :            : static GEN
     169                 :   10856748 : cycle_power_to_perm(GEN perm,GEN cy,long l)
     170                 :            : {
     171                 :   10856748 :   long lp,i,j,b, N = lg(perm), lcy = lg(cy)-1;
     172                 :            : 
     173                 :   10856748 :   lp = l % lcy;
     174         [ +  + ]:  303979501 :   for (i=1; i<N; i++) perm[i] = i;
     175         [ +  + ]:   10856748 :   if (lp)
     176                 :            :   {
     177                 :    7237811 :     pari_sp av = avma;
     178                 :    7237811 :     GEN p1 = new_chunk(N);
     179                 :    7237811 :     b = cy[1];
     180         [ +  + ]:   21713608 :     for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
     181                 :    7237811 :     perm[b] = cy[1];
     182         [ +  + ]:  202652709 :     for (i=1; i<N; i++) p1[i] = perm[i];
     183                 :            : 
     184         [ +  + ]:   10856734 :     for (j=2; j<=lp; j++) perm_mul_i(perm,p1);
     185                 :    7237811 :     avma = av;
     186                 :            :   }
     187                 :   10856748 :   return perm;
     188                 :            : }
     189                 :            : 
     190                 :            : /* image du block system D par la permutation perm */
     191                 :            : static GEN
     192                 :    1956710 : im_block_by_perm(GEN D,GEN perm)
     193                 :            : {
     194                 :    1956710 :   long i, lb = lg(D);
     195                 :    1956710 :   GEN Dn = cgetg(lb,t_VEC);
     196         [ +  + ]:   19289935 :   for (i=1; i<lb; i++) gel(Dn,i) = vecsmallpermute(perm, gel(D,i));
     197                 :    1956710 :   return Dn;
     198                 :            : }
     199                 :            : 
     200                 :            : static void
     201                 :      45717 : append(GEN D, GEN a)
     202                 :            : {
     203                 :      45717 :   long i,l = lg(D), m = lg(a);
     204                 :      45717 :   GEN x = D + (l-1);
     205         [ +  + ]:     128534 :   for (i=1; i<m; i++) gel(x,i) = gel(a,i);
     206                 :      45717 :   setlg(D, l+m-1);
     207                 :      45717 : }
     208                 :            : 
     209                 :            : static GEN
     210                 :      10640 : print_block_system(blockdata *B, GEN Y, GEN SB)
     211                 :            : {
     212                 :      10640 :   long i, j, l, ll, lp, u, v, ns, r = lg(Y), N = B->N;
     213                 :            :   long *k, *n, **e, *t;
     214                 :      10640 :   GEN D, De, Z, cyperm, perm, VOID = cgetg(1, t_VECSMALL);
     215                 :            : 
     216         [ -  + ]:      10640 :   if (DEBUGLEVEL>5) err_printf("Y = %Ps\n",Y);
     217                 :      10640 :   n = new_chunk(N+1);
     218                 :      10640 :   D = vectrunc_init(N+1);
     219                 :      10640 :   t = new_chunk(r+1);
     220                 :      10640 :   k = new_chunk(r+1);
     221                 :      10640 :   Z = cgetg(r+1, t_VEC);
     222         [ +  + ]:      56357 :   for (ns=0,i=1; i<r; i++)
     223                 :            :   {
     224                 :      45717 :     GEN Yi = gel(Y,i);
     225                 :      45717 :     long ki = 0, si = lg(Yi)-1;
     226                 :            : 
     227         [ +  + ]:     139839 :     for (j=1; j<=si; j++) { n[j] = lg(gel(Yi,j))-1; ki += n[j]; }
     228                 :      45717 :     ki /= B->size;
     229                 :      45717 :     De = cgetg(ki+1,t_VEC);
     230         [ +  + ]:     128534 :     for (j=1; j<=ki; j++) gel(De,j) = VOID;
     231         [ +  + ]:     139839 :     for (j=1; j<=si; j++)
     232                 :            :     {
     233                 :      94122 :       GEN cy = gel(Yi,j);
     234         [ +  + ]:     378140 :       for (l=1,lp=0; l<=n[j]; l++)
     235                 :            :       {
     236         [ +  + ]:     284018 :         lp++; if (lp > ki) lp = 1;
     237                 :     284018 :         gel(De,lp) = vecsmall_append(gel(De,lp), cy[l]);
     238                 :            :       }
     239                 :            :     }
     240                 :      45717 :     append(D, De);
     241 [ +  + ][ +  + ]:      45717 :     if (si>1 && ki>1)
     242                 :            :     {
     243                 :      18333 :       GEN p1 = cgetg(si,t_VEC);
     244         [ +  + ]:      54957 :       for (j=2; j<=si; j++) p1[j-1] = Yi[j];
     245                 :      18333 :       ns++;
     246                 :      18333 :       t[ns] = si-1;
     247                 :      18333 :       k[ns] = ki-1;
     248                 :      18333 :       gel(Z,ns) = p1;
     249                 :            :     }
     250                 :            :   }
     251         [ -  + ]:      10640 :   if (DEBUGLEVEL>2) err_printf("\nns = %ld\n",ns);
     252         [ +  + ]:      10640 :   if (!ns) return test_block(B, SB, D);
     253                 :            : 
     254                 :       8512 :   setlg(Z, ns+1);
     255                 :       8512 :   e = (long**)new_chunk(ns+1);
     256         [ +  + ]:      26845 :   for (i=1; i<=ns; i++)
     257                 :            :   {
     258                 :      18333 :     e[i] = new_chunk(t[i]+1);
     259         [ +  + ]:      54957 :     for (j=1; j<=t[i]; j++) e[i][j] = 0;
     260                 :            :   }
     261                 :       8512 :   cyperm= cgetg(N+1,t_VECSMALL);
     262                 :       8512 :   perm  = cgetg(N+1,t_VECSMALL); i = ns;
     263                 :            :   do
     264                 :            :   {
     265                 :    1956710 :     pari_sp av = avma;
     266         [ +  + ]:   54781881 :     for (u=1; u<=N; u++) perm[u] = u;
     267         [ +  + ]:    7247506 :     for (u=1; u<=ns; u++)
     268         [ +  + ]:   16147544 :       for (v=1; v<=t[u]; v++)
     269                 :   10856748 :         perm_mul_i(perm, cycle_power_to_perm(cyperm, gmael(Z,u,v), e[u][v]));
     270                 :    1956710 :     SB = test_block(B, SB, im_block_by_perm(D,perm));
     271                 :    1956710 :     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         [ +  + ]:    1956710 :     if (++e[ns][t[ns]] > k[ns])
     276                 :            :     {
     277                 :     652232 :       j = t[ns]-1;
     278 [ +  + ][ +  + ]:     872165 :       while (j>=1 && e[ns][j] == k[ns]) j--;
     279 [ +  + ][ +  + ]:    1099693 :       if (j >= 1) { e[ns][j]++; for (l=j+1; l<=t[ns]; l++) e[ns][l] = 0; }
     280                 :            :       else
     281                 :            :       {
     282                 :     212373 :         i = ns-1;
     283         [ +  + ]:     237874 :         while (i>=1)
     284                 :            :         {
     285                 :     229362 :           j = t[i];
     286 [ +  + ][ +  + ]:     331303 :           while (j>=1 && e[i][j] == k[i]) j--;
     287         [ +  + ]:     229362 :           if (j<1) i--;
     288                 :            :           else
     289                 :            :           {
     290                 :     203861 :             e[i][j]++;
     291         [ +  + ]:     254821 :             for (l=j+1; l<=t[i]; l++) e[i][l] = 0;
     292         [ +  + ]:     423402 :             for (ll=i+1; ll<=ns; ll++)
     293         [ +  + ]:     658602 :               for (l=1; l<=t[ll]; l++) e[ll][l] = 0;
     294                 :     203861 :             break;
     295                 :            :           }
     296                 :            :         }
     297                 :            :       }
     298                 :            :     }
     299                 :            :   }
     300         [ +  + ]:    1956710 :   while (i > 0);
     301                 :      10640 :   return SB;
     302                 :            : }
     303                 :            : 
     304                 :            : /* ALGEBRAIC PART: test potential block systems */
     305                 :            : 
     306                 :            : static GEN
     307                 :        392 : polsimplify(GEN x)
     308                 :            : {
     309                 :        392 :   long i,lx = lg(x);
     310         [ +  + ]:       2695 :   for (i=2; i<lx; i++)
     311         [ +  + ]:       2303 :     if (typ(gel(x,i)) == t_POL) gel(x,i) = constant_term(gel(x,i));
     312                 :        392 :   return x;
     313                 :            : }
     314                 :            : 
     315                 :            : /* return 0 if |g[i]| > M[i] for some i; 1 otherwise */
     316                 :            : static long
     317                 :        392 : ok_coeffs(GEN g,GEN M)
     318                 :            : {
     319                 :        392 :   long i, lg = lg(g)-1; /* g is monic, and cst term is ok */
     320         [ +  + ]:       1575 :   for (i=3; i<lg; i++)
     321         [ +  + ]:       1225 :     if (absi_cmp(gel(g,i), gel(M,i)) > 0) return 0;
     322                 :        392 :   return 1;
     323                 :            : }
     324                 :            : 
     325                 :            : /* assume x in Fq, return Tr_{Fq/Fp}(x) as a t_INT */
     326                 :            : static GEN
     327                 :       4599 : trace(GEN x, GEN Trq, GEN p)
     328                 :            : {
     329                 :            :   long i, l;
     330                 :            :   GEN s;
     331         [ -  + ]:       4599 :   if (typ(x) == t_INT) return Fp_mul(x, gel(Trq,1), p);
     332         [ -  + ]:       4599 :   l = lg(x)-1; if (l == 1) return gen_0;
     333                 :       4599 :   x++; s = mulii(gel(x,1), gel(Trq,1));
     334         [ +  + ]:      34447 :   for (i=2; i<l; i++)
     335                 :      29848 :     s = addii(s, mulii(gel(x,i), gel(Trq,i)));
     336                 :       4599 :   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                 :       1085 : poltrace(GEN x, GEN Trq, GEN p)
     342                 :            : {
     343                 :            :   long i,l;
     344                 :            :   GEN y;
     345 [ +  - ][ -  + ]:       1085 :   if (typ(x) == t_INT || varn(x) != 0) return trace(x, Trq, p);
     346                 :       1085 :   y = cgetg_copy(x, &l); y[1] = x[1];
     347         [ +  + ]:       5684 :   for (i=2; i<l; i++) gel(y,i) = trace(gel(x,i),Trq,p);
     348                 :       1085 :   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                 :        294 : chinese_retrieve_pol(GEN DATA, primedata *S, GEN listdelta)
     358                 :            : {
     359                 :        294 :   GEN interp, bezoutC, h, p = S->p, pol = FpX_red(gel(DATA,1), p);
     360                 :            :   long i, l;
     361                 :        294 :   interp = gel(DATA,9);
     362                 :        294 :   bezoutC= gel(DATA,6);
     363                 :            : 
     364                 :        294 :   h = NULL; l = lg(interp);
     365         [ +  + ]:       1379 :   for (i=1; i<l; i++)
     366                 :            :   { /* h(firstroot[i]) = listdelta[i] */
     367                 :       1085 :     GEN t = FqX_Fq_mul(gel(interp,i), gel(listdelta,i), S->T, p);
     368                 :       1085 :     t = poltrace(t, gel(S->Trk,i), p);
     369                 :       1085 :     t = FpX_mul(t, gel(bezoutC,i), p);
     370         [ +  + ]:       1085 :     h = h? FpX_add(h,t,p): t;
     371                 :            :   }
     372                 :        294 :   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                 :        294 : embedding(GEN g, GEN DATA, primedata *S, GEN den, GEN listdelta)
     382                 :            : {
     383                 :        294 :   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                 :        294 :   T   = gel(DATA,1); rt = brent_kung_optpow(degpol(T), 4, 3);
     388                 :        294 :   maxp= gel(DATA,7);
     389                 :        294 :   gp = RgX_deriv(g); av = avma;
     390                 :        294 :   w0 = chinese_retrieve_pol(DATA, S, listdelta);
     391                 :        294 :   w0_Q = centermod(gmul(w0,den), p);
     392                 :        294 :   h0 = FpXQ_inv(FpX_FpXQ_eval(gp,w0, T,p), T,p); /* = 1/g'(w0) mod (T,p) */
     393                 :        294 :   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         [ -  + ]:       1302 :     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         [ +  + ]:       1302 :     if (wpow) a = FpX_FpXQV_eval(g,wpow, T,q);
     403                 :        294 :     else      a = FpX_FpXQ_eval(g,w0, T,q); /* first time */
     404                 :            :     /* now, a = 0 (p) */
     405                 :       1302 :     a = FpXQ_mul(ZX_neg(h0), ZX_Z_divexact(a, p), T,p);
     406                 :       1302 :     w1 = ZX_add(w0, ZX_Z_mul(a, p));
     407                 :            : 
     408                 :       1302 :     w1_Q = centermod(ZX_Z_mul(w1, remii(den,q)), q);
     409         [ +  + ]:       1302 :     if (ZX_equal(w1_Q, w0_Q))
     410                 :            :     {
     411         [ +  - ]:        210 :       GEN G = is_pm1(den)? g: RgX_rescale(g,den);
     412         [ +  - ]:        210 :       if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
     413                 :            :     }
     414         [ +  + ]:       1092 :     else if (cmpii(q,maxp) > 0)
     415                 :            :     {
     416         [ +  - ]:         84 :       GEN G = is_pm1(den)? g: RgX_rescale(g,den);
     417         [ +  + ]:         84 :       if (gequal0(RgX_RgXQ_eval(G, w1_Q, T))) break;
     418         [ -  + ]:         21 :       if (DEBUGLEVEL) err_printf("coeff too big for embedding\n");
     419                 :         21 :       return NULL;
     420                 :            :     }
     421                 :       1008 :     gerepileall(av, 5, &w1,&h0,&w1_Q,&q,&p);
     422                 :       1008 :     q2 = sqri(q);
     423                 :       1008 :     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                 :       1008 :     a = FpXQ_mul(ZX_neg(h0), FpX_FpXQV_eval(gp, FpXV_red(wpow,q),T,q), T,q);
     427                 :       1008 :     a = ZX_Z_add_shallow(a, gen_1); /* 1 - h0 g'(w1) = 0 (p) */
     428                 :       1008 :     a = FpXQ_mul(h0, ZX_Z_divexact(a, p), T,p);
     429                 :       1008 :     h0 = ZX_add(h0, ZX_Z_mul(a, p));
     430                 :       1008 :     w0 = w1; w0_Q = w1_Q; p = q; q = q2;
     431                 :       1008 :   }
     432                 :        273 :   TR = gel(DATA,5);
     433         [ +  + ]:        273 :   if (!gequal0(TR)) w1_Q = RgX_translate(w1_Q, TR);
     434                 :        294 :   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                 :         84 : get_bezout(GEN pol, GEN fk, GEN p)
     441                 :            : {
     442                 :         84 :   long i, l = lg(fk);
     443                 :         84 :   GEN A, B, d, u, v, U = cgetg(l, t_VEC);
     444         [ +  + ]:        294 :   for (i=1; i<l; i++)
     445                 :            :   {
     446                 :        210 :     A = gel(fk,i);
     447                 :        210 :     B = FpX_div(pol, A, p);
     448                 :        210 :     d = FpX_extgcd(A,B,p, &u, &v);
     449         [ -  + ]:        210 :     if (degpol(d) > 0) pari_err_COPRIME("get_bezout",A,B);
     450                 :        210 :     d = constant_term(d);
     451         [ +  + ]:        210 :     if (!gequal1(d)) v = FpX_Fp_mul(v, Fp_inv(d, p), p);
     452                 :        210 :     gel(U,i) = FpX_mul(B,v, p);
     453                 :            :   }
     454                 :         84 :   return U;
     455                 :            : }
     456                 :            : 
     457                 :            : static GEN
     458                 :         84 : init_traces(GEN ff, GEN T, GEN p)
     459                 :            : {
     460                 :         84 :   long N = degpol(T),i,j,k, r = lg(ff);
     461                 :         84 :   GEN Frob = FpX_matFrobenius(T,p);
     462                 :            :   GEN y,p1,p2,Trk,pow,pow1;
     463                 :            : 
     464                 :         84 :   k = degpol(gel(ff,r-1)); /* largest degree in modular factorization */
     465                 :         84 :   pow = cgetg(k+1, t_VEC);
     466                 :         84 :   gel(pow,1) = gen_0; /* dummy */
     467                 :         84 :   gel(pow,2) = Frob;
     468                 :         84 :   pow1= cgetg(k+1, t_VEC); /* 1st line */
     469         [ +  + ]:        679 :   for (i=3; i<=k; i++)
     470                 :        595 :     gel(pow,i) = FpM_mul(gel(pow,i-1), Frob, p);
     471                 :         84 :   gel(pow1,1) = gen_0; /* dummy */
     472         [ +  + ]:        763 :   for (i=2; i<=k; i++)
     473                 :            :   {
     474                 :        679 :     p1 = cgetg(N+1, t_VEC);
     475                 :        679 :     gel(pow1,i) = p1; p2 = gel(pow,i);
     476         [ +  + ]:      11949 :     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                 :         84 :   Trk = pow; /* re-use (destroy) pow */
     481                 :         84 :   gel(Trk,1) = vec_ei(N,1);
     482         [ +  + ]:        763 :   for (i=2; i<=k; i++)
     483                 :        679 :     gel(Trk,i) = gadd(gel(Trk,i-1), gel(pow1,i));
     484                 :         84 :   y = cgetg(r, t_VEC);
     485         [ +  + ]:        294 :   for (i=1; i<r; i++) y[i] = Trk[degpol(gel(ff,i))];
     486                 :         84 :   return y;
     487                 :            : }
     488                 :            : 
     489                 :            : static void
     490                 :         84 : init_primedata(primedata *S)
     491                 :            : {
     492                 :         84 :   long i, j, l, lff = lg(S->ff), v = fetch_var(), N = degpol(S->pol);
     493                 :         84 :   GEN T, p = S->p;
     494                 :            : 
     495         [ +  + ]:         84 :   if (S->lcm == degpol(gel(S->ff,lff-1)))
     496                 :            :   {
     497                 :         77 :     T = leafcopy(gel(S->ff,lff-1));
     498                 :         77 :     setvarn(T, v);
     499                 :            :   }
     500                 :            :   else
     501                 :          7 :     T = init_Fq(p, S->lcm, v);
     502                 :         84 :   name_var(v,"y");
     503                 :         84 :   S->T = T;
     504                 :         84 :   S->firstroot = cgetg(lff, t_VECSMALL);
     505                 :         84 :   S->interp = cgetg(lff, t_VEC);
     506                 :         84 :   S->fk = cgetg(N+1, t_VEC);
     507         [ +  + ]:        294 :   for (l=1,j=1; j<lff; j++)
     508                 :            :   { /* compute roots and fix ordering (Frobenius cycles) */
     509                 :        210 :     GEN F = gel(S->ff, j), deg1 = FpX_factorff_irred(F, T,p);
     510                 :        210 :     GEN H = gel(deg1,1), a = Fq_neg(constant_term(H), T,p);
     511                 :        210 :     GEN Q = FqX_div(F, H, T,p);
     512                 :        210 :     GEN q = Fq_inv(FqX_eval(Q, a, T,p), T,p);
     513                 :        210 :     gel(S->interp,j) = FqX_Fq_mul(Q, q, T,p); /* = 1 at a, 0 at other roots */
     514                 :        210 :     S->firstroot[j] = l;
     515         [ +  + ]:       1365 :     for (i=1; i<lg(deg1); i++,l++) gel(S->fk, l) = gel(deg1, i);
     516                 :            :   }
     517                 :         84 :   S->Trk     = init_traces(S->ff, T,p);
     518                 :         84 :   S->bezoutC = get_bezout(S->pol, S->ff, p);
     519                 :         84 : }
     520                 :            : 
     521                 :            : static void
     522                 :         84 : choose_prime(primedata *S, GEN pol, GEN dpol)
     523                 :            : {
     524                 :         84 :   long i, j, k, r, lcm, oldlcm, N = degpol(pol);
     525                 :            :   ulong p, pp;
     526                 :            :   GEN Z, ff, oldff, n, oldn;
     527                 :            :   pari_sp av;
     528                 :            :   forprime_t T;
     529                 :            : 
     530                 :         84 :   u_forprime_init(&T, (N*N) >> 2, ULONG_MAX);
     531                 :         84 :   oldlcm = 0;
     532                 :         84 :   oldff = oldn = NULL; pp = 0; /* gcc -Wall */
     533                 :         84 :   av = avma;
     534 [ +  + ][ -  + ]:        672 :   for(k = 1; k < 11 || !oldlcm; k++,avma = av)
     535                 :            :   {
     536         [ -  + ]:        623 :     if (k > 5 * N) pari_err_OVERFLOW("choose_prime [too many block systems]");
     537         [ +  + ]:        630 :     do p = u_forprime_next(&T); while (!umodiu(dpol, p));
     538                 :        623 :     ff = gel(Flx_factor(ZX_to_Flx(pol,p), p), 1);
     539                 :        623 :     r = lg(ff)-1;
     540 [ +  - ][ -  + ]:        623 :     if (r == N || r >= BIL) continue;
     541                 :            : 
     542                 :        623 :     n = cgetg(r+1, t_VECSMALL); lcm = n[1] = degpol(gel(ff,1));
     543         [ +  + ]:       2919 :     for (j=2; j<=r; j++) { n[j] = degpol(gel(ff,j)); lcm = clcm(lcm, n[j]); }
     544         [ +  + ]:        623 :     if (lcm <= oldlcm) continue; /* false when oldlcm = 0 */
     545                 :            : 
     546         [ -  + ]:        168 :     if (DEBUGLEVEL) err_printf("p = %ld,\tlcm = %ld,\torbits: %Ps\n",p,lcm,n);
     547                 :        168 :     pp = p;
     548                 :        168 :     oldn = n;
     549                 :        168 :     oldff = ff;
     550         [ +  + ]:        168 :     oldlcm = lcm; if (r == 1) break;
     551                 :        133 :     av = avma;
     552                 :            :   }
     553         [ -  + ]:         84 :   if (DEBUGLEVEL) err_printf("Chosen prime: p = %ld\n", pp);
     554                 :         84 :   FlxV_to_ZXV_inplace(oldff);
     555                 :         84 :   S->ff = oldff;
     556                 :         84 :   S->lcm= oldlcm;
     557                 :         84 :   S->p  = utoipos(pp);
     558                 :         84 :   S->pol = FpX_red(pol, S->p); init_primedata(S);
     559                 :            : 
     560                 :         84 :   n = oldn; r = lg(n); Z = cgetg(r,t_VEC);
     561         [ +  + ]:        294 :   for (k=0,i=1; i<r; i++)
     562                 :            :   {
     563                 :        210 :     GEN t = cgetg(n[i]+1, t_VECSMALL); gel(Z,i) = t;
     564         [ +  + ]:       1365 :     for (j=1; j<=n[i]; j++) t[j] = ++k;
     565                 :            :   }
     566                 :         84 :   S->Z = Z;
     567                 :         84 : }
     568                 :            : 
     569                 :            : /* maxroot t_REAL */
     570                 :            : static GEN
     571                 :        273 : bound_for_coeff(long m, GEN rr, GEN *maxroot)
     572                 :            : {
     573                 :        273 :   long i,r1, lrr=lg(rr);
     574                 :        273 :   GEN p1,b1,b2,B,M, C = matpascal(m-1);
     575                 :            : 
     576         [ +  + ]:       1337 :   for (r1=1; r1 < lrr; r1++)
     577         [ +  + ]:       1323 :     if (typ(gel(rr,r1)) != t_REAL) break;
     578                 :        273 :   r1--;
     579                 :            : 
     580                 :        273 :   rr = gabs(rr,0); *maxroot = vecmax(rr);
     581         [ +  + ]:       4067 :   for (i=1; i<lrr; i++)
     582         [ +  + ]:       3794 :     if (gcmp(gel(rr,i), gen_1) < 0) gel(rr,i) = gen_1;
     583         [ +  + ]:       1337 :   for (b1=gen_1,i=1; i<=r1; i++) b1 = gmul(b1, gel(rr,i));
     584         [ +  + ]:       3003 :   for (b2=gen_1    ; i<lrr; i++) b2 = gmul(b2, gel(rr,i));
     585                 :        273 :   B = gmul(b1, gsqr(b2)); /* Mahler measure */
     586                 :        273 :   M = cgetg(m+2, t_VEC); gel(M,1) = gel(M,2) = gen_0; /* unused */
     587         [ +  + ]:       1113 :   for (i=1; i<m; i++)
     588                 :            :   {
     589                 :        840 :     p1 = gadd(gmul(gcoeff(C, m, i+1), B),/* binom(m-1, i)   */
     590                 :        840 :               gcoeff(C, m, i));          /* binom(m-1, i-1) */
     591                 :        840 :     gel(M,i+2) = ceil_safe(p1);
     592                 :            :   }
     593                 :        273 :   return M;
     594                 :            : }
     595                 :            : 
     596                 :            : /* d = requested degree for subfield. Return DATA, valid for given pol, S and d
     597                 :            :  * If DATA != NULL, translate pol [ --> pol(X+1) ] and update DATA
     598                 :            :  * 1: polynomial pol
     599                 :            :  * 2: p^e (for Hensel lifts) such that p^e > max(M),
     600                 :            :  * 3: Hensel lift to precision p^e of DATA[4]
     601                 :            :  * 4: roots of pol in F_(p^S->lcm),
     602                 :            :  * 5: number of polynomial changes (translations)
     603                 :            :  * 6: Bezout coefficients associated to the S->ff[i]
     604                 :            :  * 7: Hadamard bound for coefficients of h(x) such that g o h = 0 mod pol.
     605                 :            :  * 8: bound M for polynomials defining subfields x PD->den
     606                 :            :  * 9: *[i] = interpolation polynomial for S->ff[i] [= 1 on the first root
     607                 :            :       S->firstroot[i], 0 on the others] */
     608                 :            : static void
     609                 :        273 : compute_data(blockdata *B)
     610                 :            : {
     611                 :            :   GEN ffL, roo, pe, p1, p2, fk, fhk, MM, maxroot, pol;
     612                 :        273 :   primedata *S = B->S;
     613                 :        273 :   GEN p = S->p, T = S->T, ff = S->ff, DATA = B->DATA;
     614                 :        273 :   long i, j, l, e, N, lff = lg(ff);
     615                 :            : 
     616         [ -  + ]:        273 :   if (DEBUGLEVEL>1) err_printf("Entering compute_data()\n\n");
     617                 :        273 :   pol = B->PD->pol; N = degpol(pol);
     618                 :        273 :   roo = B->PD->roo;
     619         [ +  + ]:        273 :   if (DATA)
     620                 :            :   {
     621                 :         56 :     GEN Xm1 = gsub(pol_x(varn(pol)), gen_1);
     622                 :         56 :     GEN TR = addis(gel(DATA,5), 1);
     623                 :         56 :     GEN mTR = negi(TR), interp, bezoutC;
     624                 :            : 
     625         [ -  + ]:         56 :     if (DEBUGLEVEL>1) err_printf("... update (translate) an existing DATA\n\n");
     626                 :            : 
     627                 :         56 :     gel(DATA,5) = TR;
     628                 :         56 :     pol = RgX_translate(gel(DATA,1), gen_m1);
     629                 :         56 :     p1 = cgetg_copy(roo, &l);
     630         [ +  + ]:        651 :     for (i=1; i<l; i++) gel(p1,i) = gadd(TR, gel(roo,i));
     631                 :         56 :     roo = p1;
     632                 :            : 
     633                 :         56 :     fk = gel(DATA,4); l = lg(fk);
     634         [ +  + ]:        651 :     for (i=1; i<l; i++) gel(fk,i) = gsub(Xm1, gel(fk,i));
     635                 :            : 
     636                 :         56 :     bezoutC = gel(DATA,6); l = lg(bezoutC);
     637                 :         56 :     interp  = gel(DATA,9);
     638         [ +  + ]:        140 :     for (i=1; i<l; i++)
     639                 :            :     {
     640         [ +  - ]:         84 :       if (degpol(gel(interp,i)) > 0) /* do not turn pol_1(0) into gen_1 */
     641                 :            :       {
     642                 :         84 :         p1 = RgX_translate(gel(interp,i), gen_m1);
     643                 :         84 :         gel(interp,i) = FpXX_red(p1, p);
     644                 :            :       }
     645         [ +  + ]:         84 :       if (degpol(gel(bezoutC,i)) > 0)
     646                 :            :       {
     647                 :         49 :         p1 = RgX_translate(gel(bezoutC,i), gen_m1);
     648                 :         49 :         gel(bezoutC,i) = FpXX_red(p1, p);
     649                 :            :       }
     650                 :            :     }
     651                 :         56 :     ff = cgetg(lff, t_VEC); /* copy, do not overwrite! */
     652         [ +  + ]:        140 :     for (i=1; i<lff; i++)
     653                 :         84 :       gel(ff,i) = FpX_red(RgX_translate(gel(S->ff,i), mTR), p);
     654                 :            :   }
     655                 :            :   else
     656                 :            :   {
     657                 :        217 :     DATA = cgetg(10,t_VEC);
     658                 :        217 :     fk = S->fk;
     659                 :        217 :     gel(DATA,5) = gen_0;
     660                 :        217 :     gel(DATA,6) = leafcopy(S->bezoutC);
     661                 :        217 :     gel(DATA,9) = leafcopy(S->interp);
     662                 :            :   }
     663                 :        273 :   gel(DATA,1) = pol;
     664                 :        273 :   MM = gmul2n(bound_for_coeff(B->d, roo, &maxroot), 1);
     665                 :        273 :   gel(DATA,8) = MM;
     666                 :        273 :   e = logint(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [for d-1 test] */
     667                 :        273 :   gel(DATA,2) = pe;
     668                 :        273 :   gel(DATA,4) = roots_from_deg1(fk);
     669                 :            : 
     670                 :            :   /* compute fhk = ZpX_liftfact(pol,fk,T,p,e,pe) in 2 steps
     671                 :            :    * 1) lift in Zp to precision p^e */
     672                 :        273 :   ffL = ZpX_liftfact(pol, ff, NULL, p, e, pe);
     673                 :        273 :   fhk = NULL;
     674         [ +  + ]:        868 :   for (l=i=1; i<lff; i++)
     675                 :            :   { /* 2) lift factorization of ff[i] in Qp[X] / T */
     676                 :        595 :     GEN F, L = gel(ffL,i);
     677                 :        595 :     long di = degpol(L);
     678                 :        595 :     F = cgetg(di+1, t_VEC);
     679         [ +  + ]:       4389 :     for (j=1; j<=di; j++) F[j] = fk[l++];
     680                 :        595 :     L = ZpX_liftfact(L, F, T, p, e, pe);
     681         [ +  + ]:        595 :     fhk = fhk? shallowconcat(fhk, L): L;
     682                 :            :   }
     683                 :        273 :   gel(DATA,3) = roots_from_deg1(fhk);
     684                 :            : 
     685                 :        273 :   p1 = mulur(N, powruhalf(stor(N-1,DEFAULTPREC), N-1));
     686                 :        273 :   p2 = powru(maxroot, B->size + N*(N-1)/2);
     687                 :        273 :   p1 = divrr(mulrr(p1,p2), gsqrt(B->PD->dis,DEFAULTPREC));
     688                 :        273 :   gel(DATA,7) = mulii(shifti(ceil_safe(p1), 1), B->PD->den);
     689                 :            : 
     690         [ -  + ]:        273 :   if (DEBUGLEVEL>1) {
     691                 :          0 :     err_printf("f = %Ps\n",DATA[1]);
     692                 :          0 :     err_printf("p = %Ps, lift to p^%ld\n", p, e);
     693                 :          0 :     err_printf("2 * Hadamard bound * ind = %Ps\n",DATA[7]);
     694                 :          0 :     err_printf("2 * M = %Ps\n",DATA[8]);
     695                 :            :   }
     696         [ +  + ]:        273 :   if (B->DATA) {
     697                 :         56 :     DATA = gclone(DATA);
     698         [ +  + ]:         56 :     if (isclone(B->DATA)) gunclone(B->DATA);
     699                 :            :   }
     700                 :        273 :   B->DATA = DATA;
     701                 :        273 : }
     702                 :            : 
     703                 :            : /* g = polynomial, h = embedding. Return [[g,h]] */
     704                 :            : static GEN
     705                 :        455 : _subfield(GEN g, GEN h) { return mkvec(mkvec2(g,h)); }
     706                 :            : 
     707                 :            : /* Return a subfield, gen_0 [ change p ] or NULL [ not a subfield ] */
     708                 :            : static GEN
     709                 :    1958894 : subfield(GEN A, blockdata *B)
     710                 :            : {
     711                 :    1958894 :   long N, i, j, d, lf, m = lg(A)-1;
     712                 :            :   GEN M, pe, pol, fhk, g, e, d_1_term, delta, listdelta, whichdelta;
     713                 :    1958894 :   GEN T = B->S->T, p = B->S->p, firstroot = B->S->firstroot;
     714                 :            : 
     715                 :    1958894 :   pol= gel(B->DATA,1); N = degpol(pol); d = N/m; /* m | N */
     716                 :    1958894 :   pe = gel(B->DATA,2);
     717                 :    1958894 :   fhk= gel(B->DATA,3);
     718                 :    1958894 :   M  = gel(B->DATA,8);
     719                 :            : 
     720                 :    1958894 :   delta = cgetg(m+1,t_VEC);
     721                 :    1958894 :   whichdelta = cgetg(N+1, t_VECSMALL);
     722                 :    1958894 :   d_1_term = gen_0;
     723         [ +  + ]:   19298909 :   for (i=1; i<=m; i++)
     724                 :            :   {
     725                 :   17340015 :     GEN Ai = gel(A,i), p1 = gel(fhk,Ai[1]);
     726         [ +  + ]:   52881199 :     for (j=2; j<=d; j++)
     727                 :   35541184 :       p1 = Fq_mul(p1, gel(fhk,Ai[j]), T, pe);
     728                 :   17340015 :     gel(delta,i) = p1;
     729         [ -  + ]:   17340015 :     if (DEBUGLEVEL>5) err_printf("delta[%ld] = %Ps\n",i,p1);
     730                 :            :     /* g = prod (X - delta[i])
     731                 :            :      * if g o h = 0 (pol), we'll have h(Ai[j]) = delta[i] for all j */
     732                 :            :     /* fk[k] belongs to block number whichdelta[k] */
     733         [ +  + ]:   70221214 :     for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
     734         [ +  - ]:   17340015 :     if (typ(p1) == t_POL) p1 = constant_term(p1);
     735                 :   17340015 :     d_1_term = addii(d_1_term, p1);
     736                 :            :   }
     737                 :    1958894 :   d_1_term = centermod(d_1_term, pe); /* Tr(g) */
     738         [ +  + ]:    1958894 :   if (absi_cmp(d_1_term, gel(M,3)) > 0) {
     739         [ -  + ]:    1958502 :     if (DEBUGLEVEL>1) err_printf("d-1 test failed\n");
     740                 :    1958502 :     return NULL;
     741                 :            :   }
     742                 :        392 :   g = FqV_roots_to_pol(delta, T, pe, 0);
     743                 :        392 :   g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
     744         [ +  + ]:        392 :   if (!ok_coeffs(g,M)) {
     745         [ -  + ]:         42 :     if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
     746         [ -  + ]:         42 :     if (DEBUGLEVEL>1) err_printf("coeff too big for pol g(x)\n");
     747                 :         42 :     return NULL;
     748                 :            :   }
     749         [ +  + ]:        350 :   if (!FpX_is_squarefree(g, p)) {
     750         [ -  + ]:         56 :     if (DEBUGLEVEL>2) err_printf("pol. found = %Ps\n",g);
     751         [ -  + ]:         56 :     if (DEBUGLEVEL>1) err_printf("changing f(x): p divides disc(g)\n");
     752                 :         56 :     compute_data(B);
     753                 :         56 :     return subfield(A, B);
     754                 :            :   }
     755                 :            : 
     756                 :        294 :   lf = lg(firstroot); listdelta = cgetg(lf, t_VEC);
     757         [ +  + ]:       1379 :   for (i=1; i<lf; i++) listdelta[i] = delta[whichdelta[firstroot[i]]];
     758         [ -  + ]:        294 :   if (DEBUGLEVEL) err_printf("candidate = %Ps\n", g);
     759                 :        294 :   e = embedding(g, B->DATA, B->S, B->PD->den, listdelta);
     760         [ +  + ]:        294 :   if (!e) return NULL;
     761         [ -  + ]:        273 :   if (DEBUGLEVEL) err_printf("... OK!\n");
     762                 :    1958894 :   return _subfield(g, e);
     763                 :            : }
     764                 :            : 
     765                 :            : /* L list of current subfields, test whether potential block D is a block,
     766                 :            :  * if so, append corresponding subfield */
     767                 :            : static GEN
     768                 :    1958838 : test_block(blockdata *B, GEN L, GEN D)
     769                 :            : {
     770                 :    1958838 :   pari_sp av = avma;
     771                 :    1958838 :   GEN sub = subfield(D, B);
     772         [ +  + ]:    1958838 :   if (sub) {
     773                 :        273 :     GEN old = L;
     774         [ +  + ]:        273 :     L = gclone( L? shallowconcat(L, sub): sub );
     775         [ +  + ]:        273 :     if (old) gunclone(old);
     776                 :            :   }
     777                 :    1958838 :   avma = av; return L;
     778                 :            : }
     779                 :            : 
     780                 :            : /* subfields of degree d */
     781                 :            : static GEN
     782                 :        217 : subfields_of_given_degree(blockdata *B)
     783                 :            : {
     784                 :        217 :   pari_sp av = avma;
     785                 :            :   GEN L;
     786                 :            : 
     787         [ -  + ]:        217 :   if (DEBUGLEVEL) err_printf("\n* Look for subfields of degree %ld\n\n", B->d);
     788                 :        217 :   B->DATA = NULL; compute_data(B);
     789                 :        217 :   L = calc_block(B, B->S->Z, cgetg(1,t_VEC), NULL);
     790         [ -  + ]:        217 :   if (DEBUGLEVEL>9)
     791         [ #  # ]:          0 :     err_printf("\nSubfields of degree %ld: %Ps\n", B->d, L? L: cgetg(1,t_VEC));
     792         [ +  + ]:        217 :   if (isclone(B->DATA)) gunclone(B->DATA);
     793                 :        217 :   avma = av; return L;
     794                 :            : }
     795                 :            : 
     796                 :            : static GEN
     797                 :         91 : fix_var(GEN x, long v)
     798                 :            : {
     799                 :         91 :   long i, l = lg(x);
     800         [ +  - ]:         91 :   if (!v) return x;
     801         [ #  # ]:          0 :   for (i=1; i<l; i++) { GEN t = gel(x,i); setvarn(t[1],v); setvarn(t[2],v); }
     802                 :         91 :   return x;
     803                 :            : }
     804                 :            : 
     805                 :            : static void
     806                 :         91 : subfields_poldata(GEN T, poldata *PD)
     807                 :            : {
     808                 :            :   GEN  nf,L,dis;
     809                 :            : 
     810                 :         91 :   T = leafcopy(get_nfpol(T, &nf));
     811                 :         91 :   PD->pol = T; setvarn(T, 0);
     812         [ +  + ]:         91 :   if (nf)
     813                 :            :   {
     814                 :          7 :     PD->den = Q_denom(nf_get_zk(nf));
     815                 :          7 :     PD->roo = nf_get_roots(nf);
     816                 :          7 :     PD->dis = mulii(absi(nf_get_disc(nf)), sqri(nf_get_index(nf)));
     817                 :            :   }
     818                 :            :   else
     819                 :            :   {
     820                 :         84 :     PD->den = initgaloisborne(T,NULL,nbits2prec(bit_accuracy(ZX_max_lg(T))), &L,NULL,&dis);
     821                 :         84 :     PD->roo = L;
     822                 :         84 :     PD->dis = absi(dis);
     823                 :            :   }
     824                 :         91 : }
     825                 :            : 
     826                 :            : static GEN
     827                 :        119 : subfieldsall(GEN nf)
     828                 :            : {
     829                 :        119 :   pari_sp av = avma;
     830                 :            :   long N, ld, i, v0;
     831                 :            :   GEN G, pol, dg, LSB, NLSB;
     832                 :            :   poldata PD;
     833                 :            :   primedata S;
     834                 :            :   blockdata B;
     835                 :            : 
     836                 :            :   /* much easier if nf is Galois (WSS) */
     837                 :        119 :   G = galoisinit(nf, NULL);
     838         [ +  + ]:        119 :   if (G != gen_0)
     839                 :            :   {
     840                 :            :     GEN L, S, p;
     841                 :            :     long l;
     842                 :            : 
     843                 :         28 :     pol = get_nfpol(nf, &nf);
     844                 :         28 :     L = lift_intern( galoissubfields(G, 0, varn(pol)) );
     845                 :         28 :     l = lg(L);
     846                 :         28 :     S = cgetg(l, t_VECSMALL);
     847         [ +  + ]:        462 :     for (i=1; i<l; i++) S[i] = lg(gmael(L,i,1));
     848                 :         28 :     p = vecsmall_indexsort(S);
     849                 :         28 :     return gerepilecopy(av,  vecpermute(L, p));
     850                 :            :   }
     851                 :            : 
     852                 :         91 :   subfields_poldata(nf, &PD);
     853                 :         91 :   pol = PD.pol;
     854                 :            : 
     855                 :         91 :   v0 = varn(pol); N = degpol(pol);
     856                 :         91 :   dg = divisorsu(N); ld = lg(dg)-1;
     857         [ -  + ]:         91 :   if (DEBUGLEVEL) err_printf("\n***** Entering subfields\n\npol = %Ps\n",pol);
     858                 :            : 
     859                 :         91 :   LSB = _subfield(pol_x(0), gen_0);
     860         [ +  + ]:         91 :   if (ld > 2)
     861                 :            :   {
     862                 :         84 :     B.PD = &PD;
     863                 :         84 :     B.S  = &S;
     864                 :         84 :     B.N  = N;
     865                 :         84 :     choose_prime(&S, PD.pol, PD.dis);
     866         [ +  + ]:        301 :     for (i=ld-1; i>1; i--)
     867                 :            :     {
     868                 :        217 :       B.size  = dg[i];
     869                 :        217 :       B.d = N / B.size;
     870                 :        217 :       NLSB = subfields_of_given_degree(&B);
     871         [ +  + ]:        217 :       if (NLSB) { LSB = concat(LSB, NLSB); gunclone(NLSB); }
     872                 :            :     }
     873                 :         84 :     (void)delete_var(); /* from choose_prime */
     874                 :            :   }
     875                 :         91 :   LSB = shallowconcat(LSB, _subfield(pol, pol_x(0)));
     876         [ -  + ]:         91 :   if (DEBUGLEVEL) err_printf("\n***** Leaving subfields\n\n");
     877                 :        119 :   return fix_var(gerepilecopy(av, LSB), v0);
     878                 :            : }
     879                 :            : 
     880                 :            : GEN
     881                 :        119 : nfsubfields(GEN nf, long d)
     882                 :            : {
     883                 :        119 :   pari_sp av = avma;
     884                 :            :   long N, v0;
     885                 :            :   GEN LSB, pol, G;
     886                 :            :   poldata PD;
     887                 :            :   primedata S;
     888                 :            :   blockdata B;
     889                 :            : 
     890         [ +  - ]:        119 :   if (!d) return subfieldsall(nf);
     891                 :            : 
     892                 :          0 :   pol = get_nfpol(nf, &nf); /* in order to treat trivial cases */
     893                 :          0 :   RgX_check_ZX(pol,"nfsubfields");
     894                 :          0 :   v0 = varn(pol); N = degpol(pol);
     895         [ #  # ]:          0 :   if (d == N) return gerepilecopy(av, _subfield(pol, pol_x(v0)));
     896         [ #  # ]:          0 :   if (d == 1) return gerepilecopy(av, _subfield(pol_x(v0), pol));
     897 [ #  # ][ #  # ]:          0 :   if (d < 1 || d > N || N % d) return cgetg(1,t_VEC);
                 [ #  # ]
     898                 :            : 
     899                 :            :   /* much easier if nf is Galois (WSS) */
     900         [ #  # ]:          0 :   G = galoisinit(nf? nf: pol, NULL);
     901         [ #  # ]:          0 :   if (G != gen_0)
     902                 :            :   { /* Bingo */
     903                 :          0 :     GEN L = galoissubgroups(G), F;
     904                 :          0 :     long k,i, l = lg(L), o = N/d;
     905                 :          0 :     F = cgetg(l, t_VEC);
     906                 :          0 :     k = 1;
     907         [ #  # ]:          0 :     for (i=1; i<l; i++)
     908                 :            :     {
     909                 :          0 :       GEN H = gel(L,i);
     910         [ #  # ]:          0 :       if (group_order(H) == o)
     911                 :          0 :         gel(F,k++) = lift_intern(galoisfixedfield(G, gel(H,1), 0, v0));
     912                 :            :     }
     913                 :          0 :     setlg(F, k);
     914                 :          0 :     return gerepilecopy(av, F);
     915                 :            :   }
     916                 :            : 
     917         [ #  # ]:          0 :   subfields_poldata(nf? nf: pol, &PD);
     918                 :            : 
     919                 :          0 :   B.PD = &PD;
     920                 :          0 :   B.S  = &S;
     921                 :          0 :   B.N  = N;
     922                 :          0 :   B.d  = d;
     923                 :          0 :   B.size = N/d;
     924                 :            : 
     925                 :          0 :   choose_prime(&S, PD.pol, PD.dis);
     926                 :          0 :   LSB = subfields_of_given_degree(&B);
     927                 :          0 :   (void)delete_var(); /* from choose_prime */
     928                 :          0 :   avma = av;
     929         [ #  # ]:          0 :   if (!LSB) return cgetg(1, t_VEC);
     930                 :          0 :   G = gcopy(LSB); gunclone(LSB);
     931                 :        119 :   return fix_var(G, v0);
     932                 :            : }

Generated by: LCOV version 1.9