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 - basemath - galconj.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16375-9f41ae0) Lines: 1453 1595 91.1 %
Date: 2014-04-19 Functions: 86 91 94.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 842 1091 77.2 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000-2003  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                 :            : #include "pari.h"
      15                 :            : #include "paripriv.h"
      16                 :            : /*************************************************************************/
      17                 :            : /**                                                                     **/
      18                 :            : /**                           GALOIS CONJUGATES                         **/
      19                 :            : /**                                                                     **/
      20                 :            : /*************************************************************************/
      21                 :            : 
      22                 :         84 : static int is2sparse(GEN x)
      23                 :            : {
      24                 :         84 :   long i, l = lg(x);
      25         [ -  + ]:         84 :   if (odd(l-3)) return 0;
      26         [ +  + ]:        259 :   for(i=3; i<l; i+=2)
      27         [ +  + ]:        217 :     if (signe(gel(x,i)))
      28                 :         42 :       return 0;
      29                 :         84 :   return 1;
      30                 :            : }
      31                 :            : 
      32                 :            : static GEN
      33                 :        455 : galoisconj1(GEN nf)
      34                 :            : {
      35                 :        455 :   GEN x = get_nfpol(nf, &nf), y, z;
      36                 :        455 :   long i, lz, v = varn(x), nbmax;
      37                 :        455 :   pari_sp av = avma;
      38                 :        455 :   RgX_check_ZX(x, "nfgaloisconj");
      39                 :        455 :   nbmax = numberofconjugates(x, 2);
      40         [ +  + ]:        455 :   if (nbmax==1)
      41                 :            :   {
      42                 :        336 :     GEN res = cgetg(2,t_COL);
      43                 :        336 :     gel(res,1) = pol_x(v);
      44                 :        336 :     return res;
      45                 :            :   }
      46 [ +  + ][ +  + ]:        119 :   if (nbmax==2 && is2sparse(x))
      47                 :            :   {
      48                 :         42 :     GEN res = cgetg(3,t_COL);
      49                 :         42 :     gel(res,1) = deg1pol_shallow(gen_m1, gen_0, v);
      50                 :         42 :     gel(res,2) = pol_x(v);
      51                 :         42 :     return res;
      52                 :            :   }
      53         [ +  + ]:         77 :   if (v == 0)
      54                 :            :   {
      55                 :         42 :     long w = fetch_var();
      56         [ +  - ]:         42 :     if (nf) y = gsubst(nf, 0, pol_x(w));
      57                 :          0 :     else { y = leafcopy(x); setvarn(y, w); }
      58                 :            :   }
      59                 :            :   else
      60                 :            :   {
      61                 :         35 :     y = x;
      62                 :         35 :     x = leafcopy(x); setvarn(x, 0);
      63                 :            :   }
      64                 :         77 :   z = nfroots(y, x); lz = lg(z);
      65                 :         77 :   y = cgetg(lz, t_COL);
      66         [ +  + ]:        504 :   for (i = 1; i < lz; i++)
      67                 :            :   {
      68                 :        427 :     GEN t = lift(gel(z,i));
      69         [ +  - ]:        427 :     if (typ(t) == t_POL) setvarn(t, v);
      70                 :        427 :     gel(y,i) = t;
      71                 :            :   }
      72         [ +  + ]:         77 :   if (v == 0) delete_var();
      73                 :        455 :   return gerepileupto(av, y);
      74                 :            : }
      75                 :            : 
      76                 :            : /* nbmax: bound for the number of conjugates */
      77                 :            : static GEN
      78                 :          0 : galoisconj2pol(GEN x, long prec)
      79                 :            : {
      80                 :            :   long i, n, v, nbauto, nbmax;
      81                 :          0 :   pari_sp av = avma;
      82                 :            :   GEN y, w, polr, p1, p2;
      83                 :          0 :   n = degpol(x);
      84         [ #  # ]:          0 :   if (n <= 0) return cgetg(1, t_VEC);
      85                 :          0 :   RgX_check_ZX(x, "galoisconj2pol");
      86                 :          0 :   nbmax = numberofconjugates(x, 2);
      87                 :          0 :   polr = QX_complex_roots(x, prec);
      88                 :          0 :   p1 = gel(polr,1);
      89                 :            :   /* accuracy in decimal digits */
      90                 :          0 :   prec = (long)prec2nbits_mul(prec, LOG10_2 * 0.75);
      91                 :          0 :   nbauto = 1;
      92                 :          0 :   w = cgetg(n + 2, t_VEC);
      93                 :          0 :   gel(w,1) = gen_1;
      94         [ #  # ]:          0 :   for (i = 2; i <= n; i++) gel(w,i) = gmul(p1, gel(w,i-1));
      95                 :          0 :   v = varn(x);
      96                 :          0 :   y = cgetg(nbmax + 1, t_COL);
      97                 :          0 :   gel(y,1) = pol_x(v);
      98 [ #  # ][ #  # ]:          0 :   for (i = 2; i <= n && nbauto < nbmax; i++)
      99                 :            :   {
     100                 :          0 :     gel(w,n+1) = gel(polr,i);
     101                 :          0 :     p1 = lindep2(w, prec);
     102         [ #  # ]:          0 :     if (signe(gel(p1,n+1)))
     103                 :            :     {
     104                 :          0 :       p1[0] = evallg(n+1) | evaltyp(t_COL);
     105                 :          0 :       p2 = gdiv(RgV_to_RgX(p1, v), negi(gel(p1,n+1)));
     106         [ #  # ]:          0 :       if (gdvd(poleval(x, p2), x))
     107                 :            :       {
     108                 :          0 :         gel(y,++nbauto) = p2;
     109         [ #  # ]:          0 :         if (DEBUGLEVEL > 1) err_printf("conjugate %ld: %Ps\n", i, y[nbauto]);
     110                 :            :       }
     111                 :            :     }
     112                 :            :   }
     113         [ #  # ]:          0 :   if (nbauto < nbmax)
     114                 :          0 :     pari_warn(warner, "conjugates list may be incomplete in nfgaloisconj");
     115                 :          0 :   setlg(y, 1 + nbauto);
     116                 :          0 :   return gerepileupto(av, gen_sort(y, (void*)&gcmp, &gen_cmp_RgX));
     117                 :            : }
     118                 :            : 
     119                 :            : static GEN
     120                 :          0 : galoisconj2(GEN nf, long prec)
     121                 :            : {
     122                 :            :   long i, n, nbauto, nbmax;
     123                 :          0 :   pari_sp av = avma;
     124                 :          0 :   GEN T = get_nfpol(nf,&nf), y, w, polr, p1, p2;
     125                 :            : 
     126         [ #  # ]:          0 :   if (!nf) return galoisconj2pol(T, prec);
     127                 :          0 :   n = degpol(T);
     128         [ #  # ]:          0 :   if (n <= 0) return cgetg(1, t_VEC);
     129                 :          0 :   nbmax = numberofconjugates(T, 2);
     130                 :          0 :   y = cgetg(nbmax + 1, t_COL);
     131                 :          0 :   gel(y,1) = pol_x(varn(T));
     132                 :            :   /* accuracy in decimal digits */
     133                 :          0 :   prec = (long)prec2nbits_mul(nf_get_prec(nf), LOG10_2 * 0.75);
     134                 :          0 :   nbauto = 1;
     135                 :          0 :   polr = nf_get_allroots(nf);
     136                 :          0 :   p2 = nf_get_M(nf);
     137                 :          0 :   w = cgetg(n + 2, t_VEC);
     138         [ #  # ]:          0 :   for (i = 1; i <= n; i++) gel(w,i) = gcoeff(p2, 1, i);
     139 [ #  # ][ #  # ]:          0 :   for (i = 2; i <= n && nbauto < nbmax; i++)
     140                 :            :   {
     141                 :          0 :     gel(w,n+1) = gel(polr,i);
     142                 :          0 :     p1 = lindep2(w, prec);
     143         [ #  # ]:          0 :     if (signe(gel(p1,n+1)))
     144                 :            :     {
     145                 :          0 :       p1[0] = evallg(n+1) | evaltyp(t_COL);
     146                 :          0 :       p2 = gdiv(coltoliftalg(nf, p1), negi(gel(p1, n+1)));
     147         [ #  # ]:          0 :       if (gdvd(poleval(T, p2), T))
     148                 :            :       {
     149                 :          0 :         gel(y,++nbauto) = p2;
     150         [ #  # ]:          0 :         if (DEBUGLEVEL > 1) err_printf("conjugate %ld: %Ps\n", i, y[nbauto]);
     151                 :            :       }
     152                 :            :     }
     153                 :            :   }
     154         [ #  # ]:          0 :   if (nbauto < nbmax)
     155                 :          0 :     pari_warn(warner, "conjugates list may be incomplete in nfgaloisconj");
     156                 :          0 :   setlg(y, 1 + nbauto);
     157                 :          0 :   return gerepileupto(av, gen_sort(y, (void*)&gcmp, &gen_cmp_RgX));
     158                 :            : }
     159                 :            : /*************************************************************************/
     160                 :            : /**                                                                     **/
     161                 :            : /**                           GALOISCONJ4                               **/
     162                 :            : /**                                                                     **/
     163                 :            : /*************************************************************************/
     164                 :            : /*DEBUGLEVEL:
     165                 :            :   1: timing
     166                 :            :   2: outline
     167                 :            :   4: complete outline
     168                 :            :   6: detail
     169                 :            :   7: memory
     170                 :            :   9: complete detail
     171                 :            : */
     172                 :            : struct galois_borne {
     173                 :            :   GEN  l;
     174                 :            :   long valsol;
     175                 :            :   long valabs;
     176                 :            :   GEN  bornesol;
     177                 :            :   GEN  ladicsol;
     178                 :            :   GEN  ladicabs;
     179                 :            : };
     180                 :            : struct galois_lift {
     181                 :            :   GEN  T;
     182                 :            :   GEN  den;
     183                 :            :   GEN  p;
     184                 :            :   GEN  L;
     185                 :            :   GEN  Lden;
     186                 :            :   long e;
     187                 :            :   GEN  Q;
     188                 :            :   GEN  TQ;
     189                 :            :   struct galois_borne *gb;
     190                 :            : };
     191                 :            : struct galois_testlift {
     192                 :            :   long n;
     193                 :            :   long f;
     194                 :            :   long g;
     195                 :            :   GEN  bezoutcoeff;
     196                 :            :   GEN  pauto;
     197                 :            :   GEN  C;
     198                 :            :   GEN  Cd;
     199                 :            : };
     200                 :            : struct galois_test { /* data for permutation test */
     201                 :            :   GEN order; /* order of tests pour galois_test_perm */
     202                 :            :   GEN borne, lborne; /* coefficient bounds */
     203                 :            :   GEN ladic;
     204                 :            :   GEN PV; /* NULL or vector of test matrices (Vmatrix) */
     205                 :            :   GEN TM; /* transpose of M */
     206                 :            :   GEN L; /* p-adic roots, known mod ladic */
     207                 :            :   GEN M; /* vandermonde inverse */
     208                 :            : };
     209                 :            : /* result of the study of Frobenius degrees */
     210                 :            : enum ga_code {ga_all_normal=1,ga_ext_2=2,ga_non_wss=4};
     211                 :            : struct galois_analysis {
     212                 :            :   long p; /* prime to be lifted */
     213                 :            :   long deg; /* degree of the lift */
     214                 :            :   long ord;
     215                 :            :   long l; /* l: prime number such that T is totally split mod l */
     216                 :            :   long p4;
     217                 :            :   enum ga_code group;
     218                 :            : };
     219                 :            : struct galois_frobenius {
     220                 :            :   long p;
     221                 :            :   long fp;
     222                 :            :   long deg;
     223                 :            :   GEN Tmod;
     224                 :            :   GEN psi;
     225                 :            : };
     226                 :            : 
     227                 :            : GEN
     228                 :       2268 : vandermondeinverseprep(GEN L)
     229                 :            : {
     230                 :       2268 :   long i, j, n = lg(L);
     231                 :            :   GEN V;
     232                 :       2268 :   V = cgetg(n, t_VEC);
     233         [ +  + ]:      25753 :   for (i = 1; i < n; i++)
     234                 :            :   {
     235                 :      23485 :     pari_sp ltop = avma;
     236                 :      23485 :     GEN W = cgetg(n-1,t_VEC);
     237                 :      23485 :     long k = 1;
     238         [ +  + ]:     707588 :     for (j = 1; j < n; j++)
     239         [ +  + ]:     684103 :       if (i != j) gel(W, k++) = gsub(gel(L,i),gel(L,j));
     240                 :      23485 :     gel(V,i) = gerepileupto(ltop,divide_conquer_prod(W,&gmul));
     241                 :            :   }
     242                 :       2268 :   return V;
     243                 :            : }
     244                 :            : 
     245                 :            : /* Compute the inverse of the van der Monde matrix of T multiplied by den */
     246                 :            : GEN
     247                 :       2184 : vandermondeinverse(GEN L, GEN T, GEN den, GEN prep)
     248                 :            : {
     249                 :       2184 :   pari_sp ltop = avma;
     250                 :       2184 :   long i, n = lg(L)-1;
     251                 :            :   GEN M, P;
     252         [ -  + ]:       2184 :   if (!prep) prep = vandermondeinverseprep(L);
     253                 :       2184 :   M = cgetg(n+1, t_MAT);
     254         [ +  + ]:      24514 :   for (i = 1; i <= n; i++)
     255                 :            :   {
     256                 :      22330 :     P = RgX_Rg_div(RgX_div_by_X_x(T, gel(L,i), NULL), gel(prep,i));
     257                 :      22330 :     gel(M,i) = RgX_to_RgV(P,n);
     258                 :            :   }
     259         [ +  + ]:       2184 :   return den? gerepileupto(ltop, gmul(den, M)): gerepilecopy(ltop, M);
     260                 :            : }
     261                 :            : 
     262                 :            : /* #r = r1 + r2 */
     263                 :            : GEN
     264                 :        519 : embed_roots(GEN ro, long r1)
     265                 :            : {
     266                 :        519 :   long r2 = lg(ro)-1-r1;
     267                 :            :   GEN L;
     268         [ +  + ]:        519 :   if (!r2) L = ro;
     269                 :            :   else
     270                 :            :   {
     271                 :        463 :     long i,j, N = r1+2*r2;
     272                 :        463 :     L = cgetg(N+1, t_VEC);
     273         [ +  + ]:        587 :     for (i = 1; i <= r1; i++) gel(L,i) = gel(ro,i);
     274         [ +  + ]:       2450 :     for (j = i; j <= N; i++)
     275                 :            :     {
     276                 :       1987 :       GEN z = gel(ro,i);
     277                 :       1987 :       gel(L,j++) = z;
     278                 :       1987 :       gel(L,j++) = mkcomplex(gel(z,1), gneg(gel(z,2)));
     279                 :            :     }
     280                 :            :   }
     281                 :        519 :   return L;
     282                 :            : }
     283                 :            : GEN
     284                 :      17206 : embed_disc(GEN z, long r1, long prec)
     285                 :            : {
     286                 :      17206 :   pari_sp av = avma;
     287                 :      17206 :   GEN t = real_1(prec);
     288                 :      17206 :   long i, j, n = lg(z)-1, r2 = n-r1;
     289         [ +  + ]:      90300 :   for (i = 1; i < r1; i++)
     290                 :            :   {
     291                 :      73094 :     GEN zi = gel(z,i);
     292         [ +  + ]:     526232 :     for (j = i+1; j <= r1; j++) t = gmul(t, gsub(zi, gel(z,j)));
     293                 :            :   }
     294         [ +  + ]:      51030 :   for (j = r1+1; j <= n; j++)
     295                 :            :   {
     296                 :      33824 :     GEN zj = gel(z,j), a = gel(zj,1), b = gel(zj,2), b2 = gsqr(b);
     297         [ +  + ]:      41048 :     for (i = 1; i <= r1; i++)
     298                 :            :     {
     299                 :       7224 :       GEN zi = gel(z,i);
     300                 :       7224 :       t = gmul(t, gadd(gsqr(gsub(zi, a)), b2));
     301                 :            :     }
     302                 :      33824 :     t = gmul(t, b);
     303                 :            :   }
     304         [ +  + ]:      17206 :   if (r2) t = gmul2n(t, r2);
     305         [ +  + ]:      17206 :   if (r2 > 1)
     306                 :            :   {
     307                 :       9240 :     GEN T = real_1(prec);
     308         [ +  + ]:      33467 :     for (i = r1+1; i < n; i++)
     309                 :            :     {
     310                 :      24227 :       GEN zi = gel(z,i), a = gel(zi,1), b = gel(zi,2);
     311         [ +  + ]:      79149 :       for (j = i+1; j <= n; j++)
     312                 :            :       {
     313                 :      54922 :         GEN zj = gel(z,j), c = gel(zj,1), d = gel(zj,2);
     314                 :      54922 :         GEN f = gsqr(gsub(a,c)), g = gsqr(gsub(b,d)), h = gsqr(gadd(b,d));
     315                 :      54922 :         T = gmul(T, gmul(gadd(f,g), gadd(f,h)));
     316                 :            :       }
     317                 :            :     }
     318                 :       9240 :     t = gmul(t, T);
     319                 :            :   }
     320                 :      17206 :   t = gsqr(t);
     321         [ +  + ]:      17206 :   if (odd(r2)) t = gneg(t);
     322                 :      17206 :   return gerepileupto(av, t);
     323                 :            : }
     324                 :            : 
     325                 :            : /* Compute bound for the coefficients of automorphisms.
     326                 :            :  * T a ZX, dn a t_INT denominator or NULL */
     327                 :            : GEN
     328                 :       2268 : initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis)
     329                 :            : {
     330                 :            :   GEN L, prep, den, nf, r;
     331                 :            :   pari_timer ti;
     332                 :            : 
     333         [ -  + ]:       2268 :   if (DEBUGLEVEL>=4) timer_start(&ti);
     334                 :       2268 :   T = get_nfpol(T, &nf);
     335         [ +  + ]:       2268 :   r = nf ? nf_get_roots(nf) : NULL;
     336 [ +  + ][ +  + ]:       2268 :   if (nf &&  precision(gel(r, 1)) >= prec)
     337                 :        519 :     L = embed_roots(r, nf_get_r1(nf));
     338                 :            :   else
     339                 :       1749 :     L = QX_complex_roots(T, prec);
     340         [ -  + ]:       2268 :   if (DEBUGLEVEL>=4) timer_printf(&ti,"roots");
     341                 :       2268 :   prep = vandermondeinverseprep(L);
     342         [ +  + ]:       2268 :   if (!dn)
     343                 :            :   {
     344                 :       1148 :     GEN dis, res = divide_conquer_prod(gabs(prep,prec), mpmul);
     345                 :       1148 :     dis = ZX_disc_all(T, expi(ceil_safe(res)));
     346                 :       1148 :     den = indexpartial(T,dis);
     347         [ +  + ]:       1148 :     if (ptdis) *ptdis = dis;
     348                 :            :   }
     349                 :            :   else
     350                 :            :   {
     351 [ +  - ][ -  + ]:       1120 :     if (typ(dn) != t_INT || signe(dn) <= 0)
     352                 :          0 :       pari_err_TYPE("initgaloisborne [incorrect denominator]", dn);
     353                 :       1120 :     den = dn;
     354                 :            :   }
     355         [ +  + ]:       2268 :   if (ptprep) *ptprep = prep;
     356                 :       2268 :   *ptL = L; return den;
     357                 :            : }
     358                 :            : 
     359                 :            : /* ||| M ||| with respect to || x ||_oo. Assume M square t_MAT */
     360                 :            : GEN
     361                 :       1617 : matrixnorm(GEN M, long prec)
     362                 :            : {
     363                 :       1617 :   long i,j, n = lg(M);
     364                 :       1617 :   GEN B = real_0(prec);
     365                 :            : 
     366         [ +  + ]:      18571 :   for (i = 1; i < n; i++)
     367                 :            :   {
     368                 :      16954 :     GEN z = gabs(gcoeff(M,i,1), prec);
     369         [ +  + ]:     550396 :     for (j = 2; j < n; j++) z = gadd(z, gabs(gcoeff(M,i,j), prec));
     370         [ +  + ]:      16954 :     if (gcmp(z, B) > 0) B = z;
     371                 :            :   }
     372                 :       1617 :   return B;
     373                 :            : }
     374                 :            : 
     375                 :            : static GEN
     376                 :       1617 : galoisborne(GEN T, GEN dn, struct galois_borne *gb, long d)
     377                 :            : {
     378                 :       1617 :   pari_sp ltop = avma, av2;
     379                 :            :   GEN borne, borneroots, borneabs;
     380                 :            :   long prec;
     381                 :            :   GEN L, M, prep, den;
     382                 :            :   pari_timer ti;
     383                 :            : 
     384                 :       1617 :   prec = nbits2prec(bit_accuracy(ZX_max_lg(T)));
     385                 :       1617 :   den = initgaloisborne(T,dn,prec, &L,&prep,NULL);
     386         [ +  + ]:       1617 :   if (!dn) den = gclone(den);
     387         [ -  + ]:       1617 :   if (DEBUGLEVEL>=4) timer_start(&ti);
     388                 :       1617 :   M = vandermondeinverse(L, RgX_gtofp(T, prec), den, prep);
     389         [ -  + ]:       1617 :   if (DEBUGLEVEL>=4) timer_printf(&ti,"vandermondeinverse");
     390                 :       1617 :   borne = matrixnorm(M, prec);
     391                 :       1617 :   borneroots = gsupnorm(L, prec); /*t_REAL*/
     392                 :       1617 :   borneabs = ceil_safe(gmul(borne,gmulsg(d, powru(borneroots, d))));
     393                 :       1617 :   borneroots = ceil_safe(gmul(borne, borneroots));
     394                 :       1617 :   av2 = avma;
     395                 :            :   /*We use d-1 test, so we must overlift to 2^BITS_IN_LONG*/
     396                 :       1617 :   gb->valsol = logint(shifti(borneroots,2+BITS_IN_LONG), gb->l,NULL);
     397                 :       1617 :   gb->valabs = logint(shifti(borneabs,2), gb->l, NULL);
     398                 :       1617 :   gb->valabs = maxss(gb->valsol, gb->valabs);
     399         [ -  + ]:       1617 :   if (DEBUGLEVEL >= 4)
     400                 :          0 :     err_printf("GaloisConj:val1=%ld val2=%ld\n", gb->valsol, gb->valabs);
     401                 :       1617 :   avma = av2;
     402                 :       1617 :   gb->bornesol = gerepileuptoint(ltop, shifti(borneroots,1));
     403         [ -  + ]:       1617 :   if (DEBUGLEVEL >= 9)
     404                 :          0 :     err_printf("GaloisConj: Bound %Ps\n",borneroots);
     405                 :       1617 :   gb->ladicsol = powiu(gb->l, gb->valsol);
     406                 :       1617 :   gb->ladicabs = powiu(gb->l, gb->valabs);
     407         [ +  + ]:       1617 :   if (!dn) { dn = icopy(den); gunclone(den); }
     408                 :       1617 :   return dn;
     409                 :            : }
     410                 :            : 
     411                 :            : static GEN
     412                 :       1547 : makeLden(GEN L,GEN den, struct galois_borne *gb)
     413                 :       1547 : { return FpC_Fp_mul(L, den, gb->ladicsol); }
     414                 :            : 
     415                 :            : /* Initialize the galois_lift structure */
     416                 :            : static void
     417                 :       1645 : initlift(GEN T, GEN den, GEN p, GEN L, GEN Lden, struct galois_borne *gb, struct galois_lift *gl)
     418                 :            : {
     419                 :       1645 :   pari_sp av = avma;
     420                 :            :   long e;
     421                 :       1645 :   gl->gb = gb;
     422                 :       1645 :   gl->T = T;
     423         [ +  + ]:       1645 :   gl->den = is_pm1(den)? gen_1: den;
     424                 :       1645 :   gl->p = p;
     425                 :       1645 :   gl->L = L;
     426                 :       1645 :   gl->Lden = Lden;
     427                 :       1645 :   e = logint(shifti(gb->bornesol, 2+BITS_IN_LONG),p,NULL);
     428                 :       1645 :   avma = av;
     429         [ -  + ]:       1645 :   if (e < 2) e = 2;
     430                 :       1645 :   gl->e = e;
     431                 :       1645 :   gl->Q = powiu(p, e);
     432                 :       1645 :   gl->TQ = FpX_red(T,gl->Q);
     433                 :       1645 : }
     434                 :            : 
     435                 :            : /* Check whether f is (with high probability) a solution and compute its
     436                 :            :  * permutation */
     437                 :            : static int
     438                 :       3738 : poltopermtest(GEN f, struct galois_lift *gl, GEN pf)
     439                 :            : {
     440                 :            :   pari_sp av;
     441                 :       3738 :   GEN fx, fp, B = gl->gb->bornesol;
     442                 :            :   long i, j, ll;
     443         [ +  + ]:      33509 :   for (i = 2; i < lg(f); i++)
     444         [ +  + ]:      31563 :     if (absi_cmp(gel(f,i),B) > 0)
     445                 :            :     {
     446         [ -  + ]:       1792 :       if (DEBUGLEVEL>=4) err_printf("GaloisConj: Solution too large.\n");
     447         [ -  + ]:       1792 :       if (DEBUGLEVEL>=8) err_printf("f=%Ps\n borne=%Ps\n",f,B);
     448                 :       1792 :       return 0;
     449                 :            :     }
     450                 :       1946 :   ll = lg(gl->L);
     451                 :       1946 :   fp = const_vecsmall(ll-1, 1); /* left on stack */
     452                 :       1946 :   av = avma;
     453         [ +  + ]:      31696 :   for (i = 1; i < ll; i++, avma = av)
     454                 :            :   {
     455                 :      29750 :     fx = FpX_eval(f, gel(gl->L,i), gl->gb->ladicsol);
     456         [ +  - ]:     587188 :     for (j = 1; j < ll; j++)
     457 [ +  + ][ +  + ]:     587188 :       if (fp[j] && equalii(fx, gel(gl->Lden,j))) { pf[i]=j; fp[j]=0; break; }
     458         [ -  + ]:      29750 :     if (j == ll) return 0;
     459                 :            :   }
     460                 :       3738 :   return 1;
     461                 :            : }
     462                 :            : 
     463                 :            : static long
     464                 :       3998 : monoratlift(GEN S, GEN q, GEN qm1old,struct galois_lift *gl, GEN frob)
     465                 :            : {
     466                 :       3998 :   GEN tlift = FpX_ratlift(S,q,qm1old,qm1old,gl->den);
     467         [ +  + ]:       3998 :   if (tlift)
     468                 :            :   {
     469                 :       1005 :     pari_sp ltop = avma;
     470         [ -  + ]:       1005 :     if(DEBUGLEVEL>=4)
     471                 :          0 :       err_printf("MonomorphismLift: trying early solution %Ps\n",tlift);
     472         [ +  + ]:       1005 :     if (gl->den != gen_1) {
     473                 :        522 :       GEN N = gl->gb->ladicsol, N2 = shifti(N,-1);
     474                 :        522 :       tlift = FpX_center(FpX_red(Q_muli_to_int(tlift, gl->den), N), N,N2);
     475                 :            :     }
     476         [ +  - ]:       1005 :     if (poltopermtest(tlift, gl, frob))
     477                 :            :     {
     478         [ -  + ]:       1005 :       if(DEBUGLEVEL>=4) err_printf("MonomorphismLift: true early solution.\n");
     479                 :       1005 :       avma = ltop; return 1;
     480                 :            :     }
     481                 :          0 :     avma = ltop;
     482         [ #  # ]:          0 :     if(DEBUGLEVEL>=4) err_printf("MonomorphismLift: false early solution.\n");
     483                 :            :   }
     484                 :       3998 :   return 0;
     485                 :            : }
     486                 :            : 
     487                 :            : static GEN
     488                 :       1778 : monomorphismratlift(GEN P, GEN S, struct galois_lift *gl, GEN frob)
     489                 :            : {
     490                 :            :   pari_sp ltop, lbot;
     491                 :       1778 :   GEN Q = gl->T, p = gl->p, qold = NULL, q = p;
     492                 :       1778 :   GEN Sr, Spow, Wr = NULL, W, Prold = NULL, Pr, Qrold = NULL, Qr;
     493                 :       1778 :   long e = gl->e, level = 1, rt;
     494                 :            :   ulong mask;
     495                 :            :   GEN *gptr[2];
     496                 :            :   pari_timer ti;
     497                 :            : 
     498         [ -  + ]:       1778 :   if (DEBUGLEVEL == 1) timer_start(&ti);
     499                 :       1778 :   rt = brent_kung_optpow(degpol(P), 4, 3);
     500                 :       1778 :   mask = quadratic_prec_mask(e);
     501                 :       1778 :   Pr = FpX_red(P,q);
     502         [ -  + ]:       1778 :   Qr = (P==Q)? Pr: FpX_red(Q, q);/*A little speed up for automorphismlift*/
     503                 :       1778 :   W = FpXQ_inv(FpX_FpXQ_eval(ZX_deriv(Pr),S, Qr,q), Qr,q);
     504                 :       1778 :   gptr[0] = &S;
     505                 :       1778 :   gptr[1] = &Wr;
     506                 :            :   for (;;)
     507                 :            :   {
     508         [ -  + ]:       7081 :     if (DEBUGLEVEL>=2) { level = (level<<1) - (mask & 1); timer_start(&ti); }
     509                 :       7081 :     q = sqri(q);
     510         [ +  + ]:       7081 :     if (mask & 1) q = diviiexact(q, p);
     511                 :       7081 :     mask >>= 1;
     512                 :       7081 :     Pr = FpX_red(P, q);
     513         [ -  + ]:       7081 :     Qr = (P==Q)? Pr: FpX_red(Q, q);/*A little speed up for automorphismlift*/
     514                 :       7081 :     ltop = avma;
     515                 :       7081 :     Sr = S;
     516                 :       7081 :     Spow = FpXQ_powers(Sr, rt, Qr, q);
     517         [ +  + ]:       7081 :     if (Wr)
     518                 :            :     {
     519                 :       5303 :       W = FpXQ_mul(Wr, FpX_FpXQV_eval(ZX_deriv(Prold),FpXV_red(Spow,qold),Qrold,qold),
     520                 :            :                    Qrold,qold);
     521                 :       5303 :       W = FpXQ_mul(Wr, Fp_FpX_sub(gen_2, W, qold), Qrold,qold);
     522                 :            :     }
     523                 :       7081 :     Wr = W;
     524                 :       7081 :     S = FpXQ_mul(Wr, FpX_FpXQV_eval(Pr, Spow, Qr, q),Qr,q);
     525                 :       7081 :     lbot = avma;
     526                 :       7081 :     S = FpX_sub(Sr, S, q);
     527         [ -  + ]:       7081 :     if (DEBUGLEVEL >= 4) timer_printf(&ti, "MonomorphismLift: lift to prec %d",level);
     528         [ +  + ]:       7081 :     if (mask == 1) break;
     529                 :       6308 :     Wr = ZX_copy(Wr);
     530                 :       6308 :     gerepilemanysp(ltop, lbot, gptr, 2);
     531 [ +  + ][ +  + ]:       6308 :     if (qold && frob && monoratlift(S,q,diviiexact(qold,p),gl,frob)) return NULL;
                 [ +  + ]
     532                 :       5303 :     qold = q; Prold = Pr; Qrold = Qr;
     533                 :       5303 :   }
     534         [ -  + ]:        773 :   if (DEBUGLEVEL == 1) timer_printf(&ti, "monomorphismlift()");
     535                 :       1778 :   return S;
     536                 :            : }
     537                 :            : 
     538                 :            : /* Let T be a polynomial in Z[X] , p a prime number, S in Fp[X]/(T) so
     539                 :            :  * that T(S)=0 [p,T]. Lift S in S_0 so that T(S_0)=0 [T,p^e]
     540                 :            :  * Unclean stack */
     541                 :            : static GEN
     542                 :       1778 : automorphismlift(GEN S, struct galois_lift *gl, GEN frob)
     543                 :            : {
     544                 :       1778 :   return monomorphismratlift(gl->T, S, gl, frob);
     545                 :            : }
     546                 :            : 
     547                 :            : /* Let P be a polynomial in Z[X] , p a prime number, S in Fp[X]/(Q) so
     548                 :            :  * that T(S)=0 [p,T]. Lift S in S_0 so that T(S_0)=0 [Q,p^e]
     549                 :            :  * Unclean stack */
     550                 :            : 
     551                 :            : GEN
     552                 :          0 : monomorphismlift(GEN P, GEN S, GEN Q, GEN p, long e)
     553                 :            : {
     554                 :            :   struct galois_lift gl;
     555                 :          0 :   gl.T = Q;
     556                 :          0 :   gl.p = p;
     557                 :          0 :   gl.e = e;
     558                 :          0 :   return monomorphismratlift(P,S,&gl,NULL);
     559                 :            : }
     560                 :            : 
     561                 :            : static GEN
     562                 :       1645 : galoisdolift(struct galois_lift *gl, GEN frob)
     563                 :            : {
     564                 :       1645 :   pari_sp av = avma;
     565                 :       1645 :   GEN Tp = FpX_red(gl->T, gl->p);
     566                 :       1645 :   GEN S = FpXQ_pow(pol_x(varn(Tp)), gl->p, Tp, gl->p);
     567                 :       1645 :   return gerepileupto(av, automorphismlift(S, gl, frob));
     568                 :            : }
     569                 :            : 
     570                 :            : static void
     571                 :        539 : inittestlift(GEN plift, GEN Tmod, struct galois_lift *gl,
     572                 :            :              struct galois_testlift *gt)
     573                 :            : {
     574                 :            :   pari_timer ti;
     575                 :        539 :   gt->n = lg(gl->L) - 1;
     576                 :        539 :   gt->g = lg(Tmod) - 1;
     577                 :        539 :   gt->f = gt->n / gt->g;
     578                 :        539 :   gt->bezoutcoeff = bezout_lift_fact(gl->T, Tmod, gl->p, gl->e);
     579         [ -  + ]:        539 :   if (DEBUGLEVEL >= 2) timer_start(&ti);
     580                 :        539 :   gt->pauto = FpXQ_autpowers(plift, gt->f-1, gl->TQ, gl->Q);
     581         [ -  + ]:        539 :   if (DEBUGLEVEL >= 2) timer_printf(&ti, "Frobenius power");
     582                 :        539 : }
     583                 :            : 
     584                 :            : /* Explanation of the intheadlong technique:
     585                 :            :  * Let C be a bound, B = BITS_IN_LONG, M > C*2^B a modulus and 0 <= a_i < M for
     586                 :            :  * i=1,...,n where n < 2^B. We want to test if there exist k,l, |k| < C < M/2^B,
     587                 :            :  * such that sum a_i = k + l*M
     588                 :            :  * We write a_i*2^B/M = b_i+c_i with b_i integer and 0<=c_i<1, so that
     589                 :            :  *   sum b_i - l*2^B = k*2^B/M - sum c_i
     590                 :            :  * Since -1 < k*2^B/M < 1 and 0<=c_i<1, it follows that
     591                 :            :  *   -n-1 < sum b_i - l*2^B < 1  i.e.  -n <= sum b_i -l*2^B <= 0
     592                 :            :  * So we compute z = - sum b_i [mod 2^B] and check if 0 <= z <= n. */
     593                 :            : 
     594                 :            : /* Assume 0 <= x < mod. */
     595                 :            : long
     596                 :     960344 : intheadlong(GEN x, GEN mod)
     597                 :            : {
     598                 :     960344 :   pari_sp av = avma;
     599                 :     960344 :   long res = (long) itou(divii(shifti(x,BITS_IN_LONG),mod));
     600                 :     960344 :   avma = av; return res;
     601                 :            : }
     602                 :            : static GEN
     603                 :      21287 : vecheadlong(GEN W, GEN mod)
     604                 :            : {
     605                 :      21287 :   long i, l = lg(W);
     606                 :      21287 :   GEN V = cgetg(l, t_VECSMALL);
     607         [ +  + ]:     936292 :   for(i=1; i<l; i++) V[i] = intheadlong(gel(W,i), mod);
     608                 :      21287 :   return V;
     609                 :            : }
     610                 :            : GEN
     611                 :        945 : matheadlong(GEN W, GEN mod)
     612                 :            : {
     613                 :        945 :   long i, l = lg(W);
     614                 :        945 :   GEN V = cgetg(l,t_MAT);
     615         [ +  + ]:      22232 :   for(i=1; i<l; i++) gel(V,i) = vecheadlong(gel(W,i), mod);
     616                 :        945 :   return V;
     617                 :            : }
     618                 :            : long
     619                 :      45339 : polheadlong(GEN P, long n, GEN mod)
     620                 :            : {
     621         [ +  - ]:      45339 :   return (lg(P)>n+2)? intheadlong(gel(P,n+2),mod): 0;
     622                 :            : }
     623                 :            : 
     624                 :            : #define headlongisint(Z,N) (-(ulong)(Z)<=(ulong)(N))
     625                 :            : 
     626                 :            : static long
     627                 :        588 : frobeniusliftall(GEN sg, long el, GEN *psi, struct galois_lift *gl,
     628                 :            :                  struct galois_testlift *gt, GEN frob)
     629                 :            : {
     630                 :        588 :   pari_sp av, ltop2, ltop = avma;
     631                 :        588 :   long i,j,k, c = lg(sg)-1, n = lg(gl->L)-1, m = gt->g, d = m / c;
     632                 :            :   GEN pf, u, v, C, Cd, SG, cache;
     633                 :        588 :   long N1, N2, R1, Ni, Z, ord = gt->f, c_idx = gt->g-1;
     634                 :        588 :   long hop = 0;
     635                 :            :   GEN NN, NQ;
     636                 :            :   pari_timer ti;
     637                 :            : 
     638                 :        588 :   *psi = pf = cgetg(m, t_VECSMALL);
     639                 :        588 :   ltop2 = avma;
     640                 :        588 :   NN = diviiexact(mpfact(m), mului(c, powiu(mpfact(d), c)));
     641         [ -  + ]:        588 :   if (DEBUGLEVEL >= 4)
     642                 :          0 :     err_printf("GaloisConj:I will try %Ps permutations\n", NN);
     643                 :        588 :   N1=10000000;
     644                 :        588 :   NQ=divis_rem(NN,N1,&R1);
     645         [ -  + ]:        588 :   if (cmpiu(NQ,1000000000)>0)
     646                 :            :   {
     647                 :          0 :     pari_warn(warner,"Combinatorics too hard : would need %Ps tests!\n"
     648                 :            :         "I will skip it, but it may induce an infinite loop",NN);
     649                 :          0 :     avma = ltop; *psi = NULL; return 0;
     650                 :            :   }
     651         [ +  - ]:        588 :   N2=itos(NQ); if(!N2) N1=R1;
     652         [ -  + ]:        588 :   if (DEBUGLEVEL>=4) timer_start(&ti);
     653                 :        588 :   avma = ltop2;
     654                 :        588 :   C = gt->C;
     655                 :        588 :   Cd= gt->Cd;
     656                 :        588 :   v = FpXQ_mul(gel(gt->pauto, 1+el%ord), gel(gt->bezoutcoeff, m),gl->TQ,gl->Q);
     657         [ +  - ]:        588 :   if (gl->den != gen_1) v = FpX_Fp_mul(v, gl->den, gl->Q);
     658                 :        588 :   SG = cgetg(lg(sg),t_VECSMALL);
     659         [ +  + ]:       1806 :   for(i=1; i<lg(SG); i++) SG[i] = (el*sg[i])%ord + 1;
     660                 :        588 :   cache = cgetg(m+1,t_VECSMALL); cache[m] = polheadlong(v,1,gl->Q);
     661                 :        588 :   Z = polheadlong(v,2,gl->Q);
     662         [ +  + ]:       2618 :   for (i = 1; i < m; i++) pf[i] = 1 + i/d;
     663                 :        588 :   av = avma;
     664                 :        588 :   for (Ni = 0, i = 0; ;i++)
     665                 :            :   {
     666         [ +  + ]:    1109234 :     for (j = c_idx ; j > 0; j--)
     667                 :            :     {
     668                 :     875504 :       long h = SG[pf[j]];
     669         [ +  + ]:     875504 :       if (!mael(C,h,j))
     670                 :            :       {
     671                 :       3500 :         pari_sp av3 = avma;
     672                 :       3500 :         GEN r = FpXQ_mul(gel(gt->pauto,h), gel(gt->bezoutcoeff,j),gl->TQ,gl->Q);
     673         [ +  - ]:       3500 :         if (gl->den != gen_1) r = FpX_Fp_mul(r, gl->den, gl->Q);
     674                 :       3500 :         gmael(C,h,j) = gclone(r);
     675                 :       3500 :         mael(Cd,h,j) = polheadlong(r,1,gl->Q);
     676                 :       3500 :         avma = av3;
     677                 :            :       }
     678                 :     875504 :       uel(cache,j) = uel(cache,j+1)+umael(Cd,h,j);
     679                 :            :     }
     680         [ +  + ]:     233730 :     if (headlongisint(uel(cache,1),n))
     681                 :            :     {
     682                 :       3220 :       long ZZ = Z;
     683         [ +  + ]:      43883 :       for (j = 1; j < m; j++) ZZ += polheadlong(gmael(C,SG[pf[j]],j),2,gl->Q);
     684         [ +  + ]:       3220 :       if (headlongisint(ZZ,n))
     685                 :            :       {
     686                 :        994 :         u = v;
     687         [ +  + ]:       5817 :         for (j = 1; j < m; j++) u = ZX_add(u, gmael(C,SG[pf[j]],j));
     688                 :        994 :         u = FpX_center(FpX_red(u, gl->Q), gl->Q, shifti(gl->Q,-1));
     689         [ +  + ]:        994 :         if (poltopermtest(u, gl, frob))
     690                 :            :         {
     691         [ -  + ]:        441 :           if (DEBUGLEVEL >= 4)
     692                 :            :           {
     693                 :          0 :             timer_printf(&ti, "");
     694                 :          0 :             err_printf("GaloisConj: %d hops on %Ps tests\n",hop,addis(mulss(Ni,N1),i));
     695                 :            :           }
     696                 :        441 :           avma = ltop2; return 1;
     697                 :            :         }
     698         [ -  + ]:        553 :         if (DEBUGLEVEL >= 4) err_printf("M");
     699                 :            :       }
     700                 :       2226 :       else hop++;
     701                 :            :     }
     702 [ -  + ][ #  # ]:     233289 :     if (DEBUGLEVEL >= 4 && i % maxss(N1/20, 1) == 0)
     703                 :          0 :       timer_printf(&ti, "GaloisConj:Testing %Ps", addis(mulss(Ni,N1),i));
     704                 :     233289 :     avma = av;
     705         [ +  + ]:     233289 :     if (i == N1-1)
     706                 :            :     {
     707         [ -  + ]:        147 :       if (Ni==N2-1) N1 = R1;
     708         [ +  - ]:        147 :       if (Ni==N2) break;
     709                 :          0 :       Ni++; i = 0;
     710         [ #  # ]:          0 :       if (DEBUGLEVEL>=4) timer_start(&ti);
     711                 :            :     }
     712 [ +  - ][ +  + ]:     640332 :     for (j = 2; j < m && pf[j-1] >= pf[j]; j++)
     713                 :            :       /*empty*/; /* to kill clang Warning */
     714 [ +  + ][ +  + ]:     395045 :     for (k = 1; k < j-k && pf[k] != pf[j-k]; k++) { lswap(pf[k], pf[j-k]); }
     715         [ +  + ]:     444836 :     for (k = j - 1; pf[k] >= pf[j]; k--)
     716                 :            :       /*empty*/;
     717                 :     233142 :     lswap(pf[j], pf[k]); c_idx = j;
     718                 :     233142 :   }
     719         [ -  + ]:        147 :   if (DEBUGLEVEL>=4) err_printf("GaloisConj: not found, %d hops \n",hop);
     720                 :        588 :   *psi = NULL; avma = ltop; return 0;
     721                 :            : }
     722                 :            : 
     723                 :            : /* Compute the test matrix for the i-th line of V. Clone. */
     724                 :            : static GEN
     725                 :        945 : Vmatrix(long i, struct galois_test *td)
     726                 :            : {
     727                 :        945 :   pari_sp av = avma;
     728                 :        945 :   GEN m = gclone( matheadlong(FpC_FpV_mul(td->L, row(td->M,i), td->ladic), td->ladic));
     729                 :        945 :   avma = av; return m;
     730                 :            : }
     731                 :            : 
     732                 :            : /* Initialize galois_test */
     733                 :            : static void
     734                 :        777 : inittest(GEN L, GEN M, GEN borne, GEN ladic, struct galois_test *td)
     735                 :            : {
     736                 :        777 :   long i, n = lg(L)-1;
     737                 :        777 :   GEN p = cgetg(n+1, t_VECSMALL);
     738         [ -  + ]:        777 :   if (DEBUGLEVEL >= 8) err_printf("GaloisConj:Init Test\n");
     739                 :        777 :   td->order = p;
     740         [ +  + ]:      13076 :   for (i = 1; i <= n-2; i++) p[i] = i+2;
     741                 :        777 :   p[n-1] = 1; p[n] = 2;
     742                 :        777 :   td->borne = borne;
     743                 :        777 :   td->lborne = subii(ladic, borne);
     744                 :        777 :   td->ladic = ladic;
     745                 :        777 :   td->L = L;
     746                 :        777 :   td->M = M;
     747                 :        777 :   td->TM = shallowtrans(M);
     748                 :        777 :   td->PV = zero_zv(n);
     749                 :        777 :   gel(td->PV, 2) = Vmatrix(2, td);
     750                 :        777 : }
     751                 :            : 
     752                 :            : /* Free clones stored inside galois_test */
     753                 :            : static void
     754                 :        777 : freetest(struct galois_test *td)
     755                 :            : {
     756                 :            :   long i;
     757         [ +  + ]:      14630 :   for (i = 1; i < lg(td->PV); i++)
     758         [ +  + ]:      13853 :     if (td->PV[i]) { gunclone(gel(td->PV,i)); td->PV[i] = 0; }
     759                 :        777 : }
     760                 :            : 
     761                 :            : /* Check if the integer P seen as a p-adic number is close to an integer less
     762                 :            :  * than td->borne in absolute value */
     763                 :            : static long
     764                 :      32396 : padicisint(GEN P, struct galois_test *td)
     765                 :            : {
     766                 :      32396 :   pari_sp ltop = avma;
     767                 :      32396 :   GEN U  = modii(P, td->ladic);
     768 [ +  + ][ +  + ]:      32396 :   long r = cmpii(U, td->borne) <= 0 || cmpii(U, td->lborne) >= 0;
     769                 :      32396 :   avma = ltop; return r;
     770                 :            : }
     771                 :            : 
     772                 :            : /* Check if the permutation pf is valid according to td.
     773                 :            :  * If not, update td to make subsequent test faster (hopefully) */
     774                 :            : static long
     775                 :      68425 : galois_test_perm(struct galois_test *td, GEN pf)
     776                 :            : {
     777                 :      68425 :   pari_sp av = avma;
     778                 :      68425 :   long i, j, n = lg(td->L)-1;
     779                 :      68425 :   GEN V, P = NULL;
     780         [ +  + ]:     101633 :   for (i = 1; i < n; i++)
     781                 :            :   {
     782                 :     100366 :     long ord = td->order[i];
     783                 :     100366 :     GEN PW = gel(td->PV, ord);
     784         [ +  + ]:     100366 :     if (PW)
     785                 :            :     {
     786                 :      67970 :       ulong Z = umael(PW,1,pf[1]);
     787         [ +  + ]:    4325776 :       for (j = 2; j <= n; j++) Z += umael(PW,j,pf[j]);
     788         [ +  + ]:      67970 :       if (!headlongisint(Z,n)) break;
     789                 :            :     } else
     790                 :            :     {
     791         [ +  + ]:      32396 :       if (!P) P = vecpermute(td->L, pf);
     792                 :      32396 :       V = FpV_dotproduct(gel(td->TM,ord), P, td->ladic);
     793         [ +  + ]:      32396 :       if (!padicisint(V, td)) {
     794                 :        168 :         gel(td->PV, ord) = Vmatrix(ord, td);
     795         [ -  + ]:        168 :         if (DEBUGLEVEL >= 4) err_printf("M");
     796                 :        168 :         break;
     797                 :            :       }
     798                 :            :     }
     799                 :            :   }
     800         [ +  + ]:      68425 :   if (i == n) { avma = av; return 1; }
     801         [ -  + ]:      67158 :   if (DEBUGLEVEL >= 4) err_printf("%d.", i);
     802         [ +  + ]:      67158 :   if (i > 1)
     803                 :            :   {
     804                 :        483 :     long z = td->order[i];
     805         [ +  + ]:       1078 :     for (j = i; j > 1; j--) td->order[j] = td->order[j-1];
     806                 :        483 :     td->order[1] = z;
     807         [ -  + ]:        483 :     if (DEBUGLEVEL >= 8) err_printf("%Ps", td->order);
     808                 :            :   }
     809                 :      68425 :   avma = av; return 0;
     810                 :            : }
     811                 :            : /*Compute a*b/c when a*b will overflow*/
     812                 :            : static long
     813                 :          0 : muldiv(long a,long b,long c)
     814                 :            : {
     815                 :          0 :   return (long)((double)a*(double)b/c);
     816                 :            : }
     817                 :            : 
     818                 :            : /* F = cycle decomposition of sigma,
     819                 :            :  * B = cycle decomposition of cl(tau).
     820                 :            :  * Check all permutations pf who can possibly correspond to tau, such that
     821                 :            :  * tau*sigma*tau^-1 = sigma^s and tau^d = sigma^t, where d = ord cl(tau)
     822                 :            :  * x: vector of choices,
     823                 :            :  * G: vector allowing linear access to elts of F.
     824                 :            :  * Choices multiple of e are not changed. */
     825                 :            : static GEN
     826                 :       1344 : testpermutation(GEN F, GEN B, GEN x, long s, long e, long cut,
     827                 :            :                 struct galois_test *td)
     828                 :            : {
     829                 :       1344 :   pari_sp av, avm = avma;
     830                 :            :   long a, b, c, d, n, p1, p2, p3, p4, p5, p6, l1, l2, N1, N2, R1;
     831                 :       1344 :   long i, j, cx, hop = 0, start = 0;
     832                 :            :   GEN pf, ar, G, W, NN, NQ;
     833                 :            :   pari_timer ti;
     834         [ -  + ]:       1344 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
     835                 :       1344 :   a = lg(F)-1; b = lg(gel(F,1))-1;
     836                 :       1344 :   c = lg(B)-1; d = lg(gel(B,1))-1;
     837                 :       1344 :   n = a*b;
     838                 :       1344 :   s = (b+s) % b;
     839                 :       1344 :   pf = cgetg(n+1, t_VECSMALL);
     840                 :       1344 :   av = avma;
     841                 :       1344 :   ar = cgetg(a+2, t_VECSMALL); ar[a+1]=0;
     842                 :       1344 :   G  = cgetg(a+1, t_VECSMALL);
     843                 :       1344 :   W  = gel(td->PV, td->order[n]);
     844         [ +  + ]:      13755 :   for (cx=1, i=1, j=1; cx <= a; cx++, i++)
     845                 :            :   {
     846                 :      12411 :     gel(G,cx) = gel(F, coeff(B,i,j));
     847         [ +  + ]:      12411 :     if (i == d) { i = 0; j++; }
     848                 :            :   }
     849                 :       1344 :   NN = divis(powuu(b, c * (d - d/e)), cut);
     850         [ -  + ]:       1344 :   if (DEBUGLEVEL>=4) err_printf("GaloisConj:I will try %Ps permutations\n", NN);
     851                 :       1344 :   N1 = 1000000;
     852                 :       1344 :   NQ = divis_rem(NN,N1,&R1);
     853         [ -  + ]:       1344 :   if (cmpiu(NQ,100000000)>0)
     854                 :            :   {
     855                 :          0 :     avma = avm;
     856                 :          0 :     pari_warn(warner,"Combinatorics too hard : would need %Ps tests!\n I'll skip it but you will get a partial result...",NN);
     857                 :          0 :     return identity_perm(n);
     858                 :            :   }
     859                 :       1344 :   N2 = itos(NQ);
     860         [ +  + ]:       1484 :   for (l2 = 0; l2 <= N2; l2++)
     861                 :            :   {
     862         [ -  + ]:       1344 :     long nbiter = (l2<N2) ? N1: R1;
     863 [ -  + ][ #  # ]:       1344 :     if (DEBUGLEVEL >= 2 && N2) err_printf("%d%% ", muldiv(l2,100,N2));
     864         [ +  + ]:    4813788 :     for (l1 = 0; l1 < nbiter; l1++)
     865                 :            :     {
     866         [ +  + ]:    4813648 :       if (start)
     867                 :            :       {
     868         [ +  - ]:   12671876 :         for (i=1, j=e; i < a;)
     869                 :            :         {
     870         [ +  + ]:    7859572 :           if ((++(x[i])) != b) break;
     871                 :    3047268 :           x[i++] = 0;
     872         [ +  + ]:    3047268 :           if (i == j) { i++; j += e; }
     873                 :            :         }
     874                 :            :       }
     875                 :       1344 :       else { start=1; i = a-1; }
     876                 :            :       /* intheadlong test: overflow in + is OK, we compute mod 2^BIL */
     877         [ +  + ]:   20233794 :       for (p1 = i+1, p5 = p1%d - 1 ; p1 >= 1; p1--, p5--) /* p5 = (p1%d) - 1 */
     878                 :            :       {
     879                 :   15420146 :         ulong V = 0;
     880         [ +  + ]:   15420146 :         if (p5 == - 1) { p5 = d - 1; p6 = p1 + 1 - d; } else p6 = p1 + 1;
     881         [ +  + ]:   15420146 :         p4 = p5 ? x[p1-1] : 0;
     882                 :   15420146 :         V = 0;
     883         [ +  + ]:   51573956 :         for (p2 = 1+p4, p3 = 1 + x[p1]; p2 <= b; p2++)
     884                 :            :         {
     885                 :   36153810 :           V += umael(W,mael(G,p6,p3),mael(G,p1,p2));
     886         [ +  + ]:   36153810 :           p3 += s; if (p3 > b) p3 -= b;
     887                 :            :         }
     888         [ +  + ]:   15420146 :         p3 = 1 + x[p1] - s; if (p3 <= 0) p3 += b;
     889         [ +  + ]:   23551654 :         for (p2 = p4; p2 >= 1; p2--)
     890                 :            :         {
     891                 :    8131508 :           V += umael(W,mael(G,p6,p3),mael(G,p1,p2));
     892         [ +  + ]:    8131508 :           p3 -= s; if (p3 <= 0) p3 += b;
     893                 :            :         }
     894                 :   15420146 :         ar[p1] = ar[p1+1] + V;
     895                 :            :       }
     896         [ +  + ]:    4813648 :       if (!headlongisint(uel(ar,1),n)) continue;
     897                 :            : 
     898                 :            :       /* intheadlong succeeds. Full computation */
     899         [ +  + ]:    2114553 :       for (p1=1, p5=d; p1 <= a; p1++, p5++)
     900                 :            :       {
     901         [ +  + ]:    2046191 :         if (p5 == d) { p5 = 0; p4 = 0; } else p4 = x[p1-1];
     902         [ +  + ]:    2046191 :         if (p5 == d-1) p6 = p1+1-d; else p6 = p1+1;
     903         [ +  + ]:    5731691 :         for (p2 = 1+p4, p3 = 1 + x[p1]; p2 <= b; p2++)
     904                 :            :         {
     905                 :    3685500 :           pf[mael(G,p1,p2)] = mael(G,p6,p3);
     906         [ +  + ]:    3685500 :           p3 += s; if (p3 > b) p3 -= b;
     907                 :            :         }
     908         [ +  + ]:    2046191 :         p3 = 1 + x[p1] - s; if (p3 <= 0) p3 += b;
     909         [ +  + ]:    2671501 :         for (p2 = p4; p2 >= 1; p2--)
     910                 :            :         {
     911                 :     625310 :           pf[mael(G,p1,p2)] = mael(G,p6,p3);
     912         [ +  + ]:     625310 :           p3 -= s; if (p3 <= 0) p3 += b;
     913                 :            :         }
     914                 :            :       }
     915         [ +  + ]:      68362 :       if (galois_test_perm(td, pf))
     916                 :            :       {
     917         [ -  + ]:       1204 :         if (DEBUGLEVEL >= 1)
     918                 :            :         {
     919                 :          0 :           GEN nb = addis(mulss(l2,N1),l1);
     920                 :          0 :           timer_printf(&ti, "testpermutation(%Ps)", nb);
     921 [ #  # ][ #  # ]:          0 :           if (DEBUGLEVEL >= 2 && hop)
     922                 :          0 :             err_printf("GaloisConj:%d hop over %Ps iterations\n", hop, nb);
     923                 :            :         }
     924                 :       1204 :         avma = av; return pf;
     925                 :            :       }
     926                 :      67158 :       hop++;
     927                 :            :     }
     928                 :            :   }
     929         [ -  + ]:        140 :   if (DEBUGLEVEL >= 1)
     930                 :            :   {
     931                 :          0 :     timer_printf(&ti, "testpermutation(%Ps)", NN);
     932 [ #  # ][ #  # ]:          0 :     if (DEBUGLEVEL >= 2 && hop)
     933                 :          0 :       err_printf("GaloisConj:%d hop over %Ps iterations\n", hop, NN);
     934                 :            :   }
     935                 :       1344 :   avma = avm; return NULL;
     936                 :            : }
     937                 :            : 
     938                 :            : /* List of subgroups of (Z/mZ)^* whose order divide o, and return the list
     939                 :            :  * of their elements, sorted by increasing order */
     940                 :            : GEN
     941                 :        658 : listznstarelts(long m, long o)
     942                 :            : {
     943                 :        658 :   pari_sp av = avma;
     944                 :            :   GEN L, zn, zns, res;
     945                 :            :   long i, phi, ind, l;
     946         [ +  + ]:        658 :   if (m == 2)
     947                 :            :   {
     948                 :         35 :     res = cgetg(2, t_VEC);
     949                 :         35 :     gel(res,1) = mkvecsmall(1);
     950                 :         35 :     return res;
     951                 :            :   }
     952                 :        623 :   zn = znstar(stoi(m));
     953                 :        623 :   phi = itos(gel(zn,1));
     954                 :        623 :   o = ugcd(o, phi); /* do we impose this on input ? */
     955                 :        623 :   zns = znstar_small(zn);
     956                 :        623 :   L = cgetg(o+1, t_VEC);
     957         [ +  + ]:       1932 :   for (i=1,ind = phi; ind; ind -= phi/o, i++) /* by *decreasing* exact index */
     958                 :       1309 :     gel(L,i) = subgrouplist(gel(zn,2), mkvec(utoipos(ind)));
     959                 :        623 :   L = shallowconcat1(L); l = lg(L);
     960         [ +  + ]:       2093 :   for (i = 1; i < l; i++) gel(L,i) = znstar_hnf_elts(zns, gel(L,i));
     961                 :        658 :   return gerepilecopy(av, L);
     962                 :            : }
     963                 :            : 
     964                 :            : /* A sympol is a symmetric polynomial
     965                 :            :  *
     966                 :            :  * Currently sympol are couple of t_VECSMALL [v,w]
     967                 :            :  * v[1]...v[k], w[1]...w[k]  represent the polynomial sum(i=1,k,v[i]*s_w[i])
     968                 :            :  * where s_i(X_1,...,X_n) = sum(j=1,n,X_j^i) */
     969                 :            : 
     970                 :            : /*Return s_e*/
     971                 :            : static GEN
     972                 :       2555 : sympol_eval_newtonsum(long e, GEN O, GEN mod)
     973                 :            : {
     974                 :       2555 :   long f = lg(O), g = lg(gel(O,1)), i, j;
     975                 :       2555 :   GEN PL = cgetg(f, t_COL);
     976         [ +  + ]:      15792 :   for(i=1; i<f; i++)
     977                 :            :   {
     978                 :      13237 :     pari_sp av = avma;
     979                 :      13237 :     GEN s = gen_0;
     980         [ +  + ]:      72086 :     for(j=1; j<g; j++) s = addii(s, Fp_powu(gmael(O,i,j), (ulong)e, mod));
     981                 :      13237 :     gel(PL,i) = gerepileuptoint(av, remii(s,mod));
     982                 :            :   }
     983                 :       2555 :   return PL;
     984                 :            : }
     985                 :            : 
     986                 :            : static GEN
     987                 :       2100 : sympol_eval(GEN v, GEN NS)
     988                 :            : {
     989                 :       2100 :   pari_sp av = avma;
     990                 :            :   long i;
     991                 :       2100 :   GEN S = gen_0;
     992         [ +  + ]:       5404 :   for (i=1; i<lg(v); i++)
     993         [ +  + ]:       3304 :     if (v[i]) S = gadd(S, gmulsg(v[i], gel(NS,i)));
     994                 :       2100 :   return gerepileupto(av, S);
     995                 :            : }
     996                 :            : 
     997                 :            : /* Let sigma be an automorphism of L (as a polynomial with rational coefs)
     998                 :            :  * Let 'sym' be a symmetric polynomial defining alpha in L.
     999                 :            :  * We have alpha = sym(x,sigma(x),,,sigma^(g-1)(x)). Compute alpha mod p */
    1000                 :            : static GEN
    1001                 :        770 : sympol_aut_evalmod(GEN sym, long g, GEN sigma, GEN Tp, GEN p)
    1002                 :            : {
    1003                 :        770 :   pari_sp ltop=avma;
    1004                 :        770 :   long i, j, npows = brent_kung_optpow(degpol(Tp)-1, g-1, 1);
    1005                 :        770 :   GEN s, f, pows, v = zv_to_ZV(gel(sym,1)), w = zv_to_ZV(gel(sym,2));
    1006                 :        770 :   sigma = RgX_to_FpX(sigma, p);
    1007                 :        770 :   pows  = FpXQ_powers(sigma,npows,Tp,p);
    1008                 :        770 :   f = pol_x(varn(sigma));
    1009                 :        770 :   s = pol_0(varn(sigma));
    1010         [ +  + ]:       4382 :   for(i=1; i<=g;i++)
    1011                 :            :   {
    1012         [ +  + ]:       3612 :     if (i > 1) f = FpX_FpXQV_eval(f,pows,Tp,p);
    1013         [ +  + ]:       7693 :     for(j=1; j<lg(v); j++)
    1014                 :       4081 :       s = FpX_add(s, FpX_Fp_mul(FpXQ_pow(f,gel(w,j),Tp,p),gel(v,j),p),p);
    1015                 :            :   }
    1016                 :        770 :   return gerepileupto(ltop, s);
    1017                 :            : }
    1018                 :            : 
    1019                 :            : /* Let Sp be as computed with sympol_aut_evalmod
    1020                 :            :  * Let Tmod be the factorisation of T mod p.
    1021                 :            :  * Return the factorisation of the minimal polynomial of S mod p w.r.t. Tmod */
    1022                 :            : static GEN
    1023                 :        770 : fixedfieldfactmod(GEN Sp, GEN p, GEN Tmod)
    1024                 :            : {
    1025                 :        770 :   long i, l = lg(Tmod);
    1026                 :        770 :   GEN F = cgetg(l,t_VEC);
    1027         [ +  + ]:       3661 :   for(i=1; i<l; i++)
    1028                 :            :   {
    1029                 :       2891 :     GEN Ti = gel(Tmod,i);
    1030                 :       2891 :     gel(F,i) = FpXQ_minpoly(FpX_rem(Sp,Ti,p), Ti,p);
    1031                 :            :   }
    1032                 :        770 :   return F;
    1033                 :            : }
    1034                 :            : 
    1035                 :            : static GEN
    1036                 :       1631 : fixedfieldsurmer(GEN mod, GEN l, GEN p, long v, GEN NS, GEN W)
    1037                 :            : {
    1038                 :       1631 :   const long step=3;
    1039                 :       1631 :   long i, j, n = lg(W)-1, m = 1L<<((n-1)<<1);
    1040                 :       1631 :   GEN sym = cgetg(n+1,t_VECSMALL), mod2 = shifti(mod,-1);
    1041         [ +  + ]:       1918 :   for (j=1;j<n;j++) sym[j] = step;
    1042                 :       1631 :   sym[n] = 0;
    1043         [ -  + ]:       1631 :   if (DEBUGLEVEL>=4) err_printf("FixedField: Weight: %Ps\n",W);
    1044         [ +  + ]:       2212 :   for (i=0; i<m; i++)
    1045                 :            :   {
    1046                 :       2100 :     pari_sp av = avma;
    1047                 :            :     GEN L, P;
    1048         [ +  + ]:       2457 :     for (j=1; sym[j]==step; j++) sym[j]=0;
    1049                 :       2100 :     sym[j]++;
    1050         [ -  + ]:       2100 :     if (DEBUGLEVEL>=6) err_printf("FixedField: Sym: %Ps\n",sym);
    1051                 :       2100 :     L = sympol_eval(sym,NS);
    1052         [ +  + ]:       2100 :     if (!vec_is1to1(FpC_red(L,l))) continue;
    1053                 :       1687 :     P = FpX_center(FpV_roots_to_pol(L,mod,v),mod,mod2);
    1054 [ +  + ][ +  + ]:       1687 :     if (!p || FpX_is_squarefree(P,p)) return mkvec3(mkvec2(sym,W),L,P);
    1055                 :        168 :     avma = av;
    1056                 :            :   }
    1057                 :       1631 :   return NULL;
    1058                 :            : }
    1059                 :            : 
    1060                 :            : /*Check whether the line of NS are pair-wise distinct.*/
    1061                 :            : static long
    1062                 :       1757 : sympol_is1to1_lg(GEN NS, long n)
    1063                 :            : {
    1064                 :       1757 :   long i, j, k, l = lgcols(NS);
    1065         [ +  + ]:      10871 :   for (i=1; i<l; i++)
    1066         [ +  + ]:      59731 :     for(j=i+1; j<l; j++)
    1067                 :            :     {
    1068         [ +  + ]:      52059 :       for(k=1; k<n; k++)
    1069         [ +  + ]:      51933 :         if (!equalii(gmael(NS,k,j),gmael(NS,k,i))) break;
    1070         [ +  + ]:      50617 :       if (k>=n) return 0;
    1071                 :            :     }
    1072                 :       1757 :   return 1;
    1073                 :            : }
    1074                 :            : 
    1075                 :            : /* Let O a set of orbits of roots (see fixedfieldorbits) modulo mod,
    1076                 :            :  * l | mod and p two prime numbers. Return a vector [sym,s,P] where:
    1077                 :            :  * sym is a sympol, s is the set of images of sym on O and
    1078                 :            :  * P is the polynomial with roots s. */
    1079                 :            : static GEN
    1080                 :       1519 : fixedfieldsympol(GEN O, GEN mod, GEN l, GEN p, long v)
    1081                 :            : {
    1082                 :       1519 :   pari_sp ltop=avma;
    1083                 :       1519 :   const long n=(BITS_IN_LONG>>1)-1;
    1084                 :       1519 :   GEN NS = cgetg(n+1,t_MAT), sym = NULL, W = cgetg(n+1,t_VECSMALL);
    1085                 :       1519 :   long i, e=1;
    1086         [ -  + ]:       1519 :   if (DEBUGLEVEL>=4)
    1087                 :          0 :     err_printf("FixedField: Size: %ldx%ld\n",lg(O)-1,lg(gel(O,1))-1);
    1088 [ +  + ][ +  - ]:       3276 :   for (i=1; !sym && i<=n; i++)
    1089                 :            :   {
    1090                 :       1757 :     GEN L = sympol_eval_newtonsum(e++, O, mod);
    1091         [ +  + ]:       1757 :     if (lg(O)>2)
    1092         [ +  + ]:       2513 :       while (vec_isconst(L)) L = sympol_eval_newtonsum(e++, O, mod);
    1093                 :       1757 :     W[i] = e-1; gel(NS,i) = L;
    1094         [ +  + ]:       1757 :     if (sympol_is1to1_lg(NS,i+1))
    1095                 :       1631 :       sym = fixedfieldsurmer(mod,l,p,v,NS,vecsmall_shorten(W,i));
    1096                 :            :   }
    1097         [ -  + ]:       1519 :   if (!sym) pari_err_BUG("fixedfieldsympol [p too small]");
    1098         [ -  + ]:       1519 :   if (DEBUGLEVEL>=2) err_printf("FixedField: Found: %Ps\n",gel(sym,1));
    1099                 :       1519 :   return gerepilecopy(ltop,sym);
    1100                 :            : }
    1101                 :            : 
    1102                 :            : /* Let O a set of orbits as indices and L the corresponding roots.
    1103                 :            :  * Return the set of orbits as roots. */
    1104                 :            : static GEN
    1105                 :       1519 : fixedfieldorbits(GEN O, GEN L)
    1106                 :            : {
    1107                 :       1519 :   GEN S = cgetg(lg(O), t_MAT);
    1108                 :            :   long i;
    1109         [ +  + ]:       9233 :   for (i = 1; i < lg(O); i++) gel(S,i) = vecpermute(L, gel(O,i));
    1110                 :       1519 :   return S;
    1111                 :            : }
    1112                 :            : 
    1113                 :            : static GEN
    1114                 :        539 : fixedfieldinclusion(GEN O, GEN PL)
    1115                 :            : {
    1116                 :        539 :   long i, j, f = lg(O)-1, g = lg(gel(O,1))-1;
    1117                 :        539 :   GEN S = cgetg(f*g + 1, t_COL);
    1118         [ +  + ]:       3304 :   for (i = 1; i <= f; i++)
    1119                 :            :   {
    1120                 :       2765 :     GEN Oi = gel(O,i);
    1121         [ +  + ]:      10549 :     for (j = 1; j <= g; j++) gel(S, Oi[j]) = gel(PL, i);
    1122                 :            :   }
    1123                 :        539 :   return S;
    1124                 :            : }
    1125                 :            : 
    1126                 :            : /*Usually mod > den so there is no need to reduce it.*/
    1127                 :            : GEN
    1128                 :       1659 : vandermondeinversemod(GEN L, GEN T, GEN den, GEN mod)
    1129                 :            : {
    1130                 :            :   pari_sp av;
    1131                 :       1659 :   long i, n = lg(L);
    1132                 :       1659 :   GEN P, Tp, M = cgetg(n, t_MAT);
    1133                 :       1659 :   av = avma;
    1134                 :       1659 :   Tp = gclone(FpX_deriv(T,mod)); /*clone*/
    1135                 :       1659 :   avma = av;
    1136         [ +  + ]:      19453 :   for (i = 1; i < n; i++)
    1137                 :            :   {
    1138                 :            :     GEN z;
    1139                 :      17794 :     av = avma;
    1140                 :      17794 :     z = Fp_inv(FpX_eval(Tp, gel(L,i),mod),mod);
    1141                 :      17794 :     z = Fp_mul(den,z,mod);
    1142                 :      17794 :     P = FpX_Fp_mul(FpX_div_by_X_x(T, gel(L,i), mod, NULL), z, mod);
    1143                 :      17794 :     gel(M,i) = gerepilecopy(av, RgX_to_RgV(P, n-1));
    1144                 :            :   }
    1145                 :       1659 :   gunclone(Tp); /*unclone*/
    1146                 :       1659 :   return M;
    1147                 :            : }
    1148                 :            : 
    1149                 :            : /* Polynomial associated to a vector of conjugates. Not stack clean */
    1150                 :            : static GEN
    1151                 :       5656 : vectopol(GEN v, GEN M, GEN den , GEN mod, GEN mod2, long x)
    1152                 :            : {
    1153                 :       5656 :   long l = lg(v)+1, i;
    1154                 :       5656 :   GEN z = cgetg(l,t_POL);
    1155                 :       5656 :   z[1] = evalsigne(1)|evalvarn(x);
    1156         [ +  + ]:      62706 :   for (i=2; i<l; i++)
    1157                 :      57050 :     gel(z,i) = gdiv(centermodii(ZMrow_ZC_mul(M,v,i-1), mod, mod2), den);
    1158                 :       5656 :   return normalizepol_lg(z, l);
    1159                 :            : }
    1160                 :            : 
    1161                 :            : /* Polynomial associate to a permutation of the roots. Not stack clean */
    1162                 :            : static GEN
    1163                 :       4893 : permtopol(GEN p, GEN L, GEN M, GEN den, GEN mod, GEN mod2, long x)
    1164                 :            : {
    1165         [ -  + ]:       4893 :   if (lg(p) != lg(L)) pari_err_TYPE("permtopol [permutation]", p);
    1166                 :       4893 :   return vectopol(vecpermute(L,p), M, den, mod, mod2, x);
    1167                 :            : }
    1168                 :            : 
    1169                 :            : static GEN
    1170                 :        847 : galoisgrouptopol(GEN res, GEN L, GEN M, GEN den, GEN mod, long v)
    1171                 :            : {
    1172                 :        847 :   long i, l = lg(res);
    1173                 :        847 :   GEN mod2 = shifti(mod,-1), aut = cgetg(l, t_COL);
    1174         [ +  + ]:       3465 :   for (i = 1; i < l; i++)
    1175                 :            :   {
    1176         [ -  + ]:       2618 :     if (DEBUGLEVEL>=6) err_printf("%d ",i);
    1177                 :       2618 :     gel(aut,i) = permtopol(gel(res,i), L, M, den, mod, mod2, v);
    1178                 :            :   }
    1179                 :        847 :   return aut;
    1180                 :            : }
    1181                 :            : 
    1182                 :            : static void
    1183                 :        546 : notgalois(long p, struct galois_analysis *ga)
    1184                 :            : {
    1185         [ -  + ]:        546 :   if (DEBUGLEVEL >= 2) err_printf("GaloisAnalysis:non Galois for p=%ld\n", p);
    1186                 :        546 :   ga->p = p;
    1187                 :        546 :   ga->deg = 0;
    1188                 :        546 : }
    1189                 :            : 
    1190                 :            : #define numberof(x) (long)(sizeof(x) / sizeof((x)[0]))
    1191                 :            : 
    1192                 :            : /*Gather information about the group*/
    1193                 :            : static long
    1194                 :       2114 : init_group(long n, long np, GEN Fp, GEN Fe, long *porder)
    1195                 :            : {
    1196                 :            :   /*TODO: complete the table to at least 200*/
    1197                 :       2114 :   const long prim_nonss_orders[] = { 36,48,56,60,72,75,80,96,108,120,132 };
    1198                 :       2114 :   long i, phi_order = 1, order = 1, group = 0;
    1199                 :            : 
    1200                 :            :  /* non-WSS groups of this order? */
    1201         [ +  + ]:      25228 :   for (i=0; i < numberof(prim_nonss_orders); i++)
    1202         [ +  + ]:      23128 :     if (n % prim_nonss_orders[i] == 0) { group |= ga_non_wss; break; }
    1203 [ +  + ][ +  + ]:       2114 :   if (np == 2 && Fp[2] == 3 && Fe[2] == 1 && Fe[1] > 2) group |= ga_ext_2;
         [ +  + ][ +  + ]
    1204                 :            : 
    1205         [ +  + ]:       3430 :   for (i = np; i > 0; i--)
    1206                 :            :   {
    1207                 :       2457 :     long p = Fp[i];
    1208         [ +  + ]:       2457 :     if (phi_order % p == 0) { group |= ga_all_normal; break; }
    1209                 :       2142 :     order *= p; phi_order *= p-1;
    1210         [ +  + ]:       2142 :     if (Fe[i] > 1) break;
    1211                 :            :   }
    1212                 :       2114 :   *porder = order; return group;
    1213                 :            : }
    1214                 :            : 
    1215                 :            : /*is a "better" than b ? (if so, update karma) */
    1216                 :            : static int
    1217                 :      11935 : improves(long a, long b, long plift, long p, long n, long *karma)
    1218                 :            : {
    1219 [ +  + ][ +  + ]:      11935 :   if (!plift || a > b) { *karma = ugcd(p-1,n); return 1; }
    1220         [ +  + ]:      11375 :   if (a == b) {
    1221                 :       9632 :     long k = ugcd(p-1,n);
    1222         [ +  + ]:       9632 :     if (k > *karma) { *karma = k; return 1; }
    1223                 :            :   }
    1224                 :      11935 :   return 0; /* worse */
    1225                 :            : }
    1226                 :            : 
    1227                 :            : /* return 0 if not galois or not wss */
    1228                 :            : static int
    1229                 :       2114 : galoisanalysis(GEN T, struct galois_analysis *ga, long calcul_l)
    1230                 :            : {
    1231                 :       2114 :   pari_sp ltop = avma, av;
    1232                 :       2114 :   long group, linf, n, p, i, karma = 0;
    1233                 :            :   GEN F, Fp, Fe, Fpe, O;
    1234                 :            :   long np, order, plift, nbmax, nbtest, deg;
    1235                 :            :   pari_timer ti;
    1236                 :            :   forprime_t S;
    1237         [ -  + ]:       2114 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
    1238                 :       2114 :   n = degpol(T);
    1239                 :       2114 :   O = zero_zv(n);
    1240                 :       2114 :   F = factoru_pow(n);
    1241                 :       2114 :   Fp = gel(F,1); np = lg(Fp)-1;
    1242                 :       2114 :   Fe = gel(F,2);
    1243                 :       2114 :   Fpe= gel(F,3);
    1244                 :       2114 :   group = init_group(n, np, Fp, Fe, &order);
    1245                 :            : 
    1246                 :            :   /*Now we study the orders of the Frobenius elements*/
    1247                 :       2114 :   deg = Fp[np]; /* largest prime | n */
    1248                 :       2114 :   plift = 0;
    1249                 :       2114 :   nbtest = 0;
    1250                 :       2114 :   nbmax = 8+(n>>1);
    1251                 :       2114 :   u_forprime_init(&S, n*maxss(expu(n)-3, 2), ULONG_MAX);
    1252                 :       2114 :   av = avma;
    1253 [ +  + ][ +  + ]:      23030 :   while (!plift || (nbtest < nbmax && (nbtest <=8 || order < (n>>1)))
         [ +  + ][ +  + ]
    1254 [ +  + ][ +  + ]:       1568 :                 || (n == 24 && O[6] == 0 && O[4] == 0)
                 [ -  + ]
    1255 [ +  + ][ -  + ]:       1568 :                 || ((group&ga_non_wss) && order == Fp[np]))
    1256                 :            :   {
    1257                 :      21462 :     long d, o, norm_o = 1;
    1258                 :            :     GEN D, Tp;
    1259                 :            : 
    1260 [ +  + ][ +  - ]:      21462 :     if ((group&ga_non_wss) && nbtest >= 3*nbmax) break; /* in all cases */
    1261                 :      21462 :     nbtest++; avma = av;
    1262                 :      21462 :     p = u_forprime_next(&S);
    1263         [ -  + ]:      21462 :     if (!p) pari_err_OVERFLOW("galoisanalysis [ran out of primes]");
    1264                 :      21462 :     Tp = ZX_to_Flx(T,p);
    1265 [ +  + ][ +  + ]:      21462 :     if (!Flx_is_squarefree(Tp,p)) { if (!--nbtest) nbtest = 1; continue; }
    1266                 :            : 
    1267                 :      19901 :     D = Flx_nbfact_by_degree(Tp, &d, p);
    1268                 :      19901 :     o = n / d; /* d factors, all should have degree o */
    1269         [ +  + ]:      19901 :     if (D[o] != d) { notgalois(p, ga); avma = ltop; return 0; }
    1270                 :            : 
    1271         [ +  + ]:      19355 :     if (!O[o]) O[o] = p;
    1272         [ +  + ]:      19355 :     if (o % deg) goto ga_end; /* NB: deg > 1 */
    1273 [ +  + ][ +  + ]:      13531 :     if ((group&ga_all_normal) && o < order) goto ga_end;
    1274                 :            : 
    1275                 :            :     /*Frob_p has order o > 1, find a power which generates a normal subgroup*/
    1276         [ +  + ]:      13426 :     if (o * Fp[1] >= n)
    1277                 :       7105 :       norm_o = o; /*subgroups of smallest index are normal*/
    1278                 :            :     else
    1279                 :            :     {
    1280         [ +  - ]:       6895 :       for (i = np; i > 0; i--)
    1281                 :            :       {
    1282         [ +  + ]:       6895 :         if (o % Fpe[i]) break;
    1283                 :        574 :         norm_o *= Fpe[i];
    1284                 :            :       }
    1285                 :            :     }
    1286                 :            :     /* Frob_p^(o/norm_o) generates a normal subgroup of order norm_o */
    1287         [ +  + ]:      13426 :     if (norm_o != 1)
    1288                 :            :     {
    1289 [ +  + ][ +  + ]:       7679 :       if (!(group&ga_all_normal) || o > order)
    1290                 :       1491 :         karma = ugcd(p-1,n);
    1291         [ +  + ]:       6188 :       else if (!improves(norm_o, deg, plift,p,n, &karma)) goto ga_end;
    1292                 :            :       /* karma0=0, deg0<=norm_o -> the first improves() returns 1 */
    1293                 :       2163 :       deg = norm_o; group |= ga_all_normal; /* STORE */
    1294                 :            :     }
    1295         [ -  + ]:       5747 :     else if (group&ga_all_normal) goto ga_end;
    1296         [ +  + ]:       5747 :     else if (!improves(o, order, plift,p,n, &karma)) goto ga_end;
    1297                 :            : 
    1298                 :       2751 :     order = o; plift = p; /* STORE */
    1299                 :            :     ga_end:
    1300         [ -  + ]:      19355 :     if (DEBUGLEVEL >= 5)
    1301                 :      19355 :       err_printf("GaloisAnalysis:Nbtest=%ld,p=%ld,o=%ld,n_o=%d,best p=%ld,ord=%ld,k=%ld\n", nbtest, p, o, norm_o, plift, order,karma);
    1302                 :            :   }
    1303                 :            :   /* To avoid looping on non-wss group.
    1304                 :            :    * TODO: check for large groups. Would it be better to disable this check if
    1305                 :            :    * we are in a good case (ga_all_normal && !(ga_ext_2) (e.g. 60)) ?*/
    1306                 :       1568 :   ga->p = plift;
    1307 [ +  - ][ +  + ]:       1568 :   if (!plift || ((group&ga_non_wss) && order == Fp[np]))
                 [ -  + ]
    1308                 :            :   {
    1309                 :          0 :     pari_warn(warner,"Galois group almost certainly not weakly super solvable");
    1310                 :          0 :     return 0;
    1311                 :            :   }
    1312                 :            :   /*linf=(n*(n-1))>>2;*/
    1313                 :       1568 :   linf = n;
    1314 [ +  + ][ +  + ]:       1568 :   if (calcul_l && O[1] <= linf)
    1315                 :            :   {
    1316                 :            :     pari_sp av2;
    1317                 :            :     forprime_t S2;
    1318                 :            :     ulong p;
    1319                 :        539 :     u_forprime_init(&S2, linf+1,ULONG_MAX);
    1320                 :        539 :     av2 = avma;
    1321         [ +  - ]:      41923 :     while ((p = u_forprime_next(&S2)))
    1322                 :            :     { /*find a totally split prime l > linf*/
    1323                 :      41923 :       GEN Tp = ZX_to_Flx(T, p);
    1324                 :      41923 :       long nb = Flx_nbroots(Tp, p);
    1325         [ +  + ]:      41923 :       if (nb == n) { O[1] = p; break; }
    1326 [ +  + ][ -  + ]:      41384 :       if (nb && Flx_is_squarefree(Tp,p)) {
    1327                 :          0 :         notgalois(p,ga);
    1328                 :          0 :         avma = ltop; return 0;
    1329                 :            :       }
    1330                 :      41384 :       avma = av2;
    1331                 :            :     }
    1332         [ -  + ]:        539 :     if (!p) pari_err_OVERFLOW("galoisanalysis [ran out of primes]");
    1333                 :            :   }
    1334                 :       1568 :   ga->group = (enum ga_code)group;
    1335                 :       1568 :   ga->deg = deg;
    1336                 :       1568 :   ga->ord = order;
    1337                 :       1568 :   ga->l  = O[1];
    1338         [ +  + ]:       1568 :   ga->p4 = n >= 4 ? O[4] : 0;
    1339         [ -  + ]:       1568 :   if (DEBUGLEVEL >= 4)
    1340                 :          0 :     err_printf("GaloisAnalysis:p=%ld l=%ld group=%ld deg=%ld ord=%ld\n",
    1341                 :          0 :                plift, O[1], group, deg, order);
    1342         [ -  + ]:       1568 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "galoisanalysis()");
    1343                 :       2114 :   avma = ltop; return 1;
    1344                 :            : }
    1345                 :            : 
    1346                 :            : static GEN
    1347                 :         21 : a4galoisgen(struct galois_test *td)
    1348                 :            : {
    1349                 :         21 :   const long n = 12;
    1350                 :         21 :   pari_sp ltop = avma, av, av2;
    1351                 :         21 :   long i, j, k, N, hop = 0;
    1352                 :            :   GEN MT, O,O1,O2,O3, ar, mt, t, u, res, orb, pft, pfu, pfv;
    1353                 :            : 
    1354                 :         21 :   res = cgetg(3, t_VEC);
    1355                 :         21 :   pft = cgetg(n+1, t_VECSMALL);
    1356                 :         21 :   pfu = cgetg(n+1, t_VECSMALL);
    1357                 :         21 :   pfv = cgetg(n+1, t_VECSMALL);
    1358                 :         21 :   gel(res,1) = mkvec3(pft,pfu,pfv);
    1359                 :         21 :   gel(res,2) = mkvecsmall3(2,2,3);
    1360                 :         21 :   av = avma;
    1361                 :         21 :   ar = cgetg(5, t_VECSMALL);
    1362                 :         21 :   mt = gel(td->PV, td->order[n]);
    1363                 :         21 :   t = identity_perm(n) + 1; /* Sorry for this hack */
    1364                 :         21 :   u = cgetg(n+1, t_VECSMALL) + 1; /* too lazy to correct */
    1365                 :         21 :   MT = cgetg(n+1, t_MAT);
    1366         [ +  + ]:        273 :   for (j = 1; j <= n; j++) gel(MT,j) = cgetg(n+1, t_VECSMALL);
    1367         [ +  + ]:        273 :   for (j = 1; j <= n; j++)
    1368         [ +  + ]:       1638 :     for (i = 1; i < j; i++)
    1369                 :       1386 :       ucoeff(MT,i,j) = ucoeff(MT,j,i) = ucoeff(mt,i,j)+ucoeff(mt,j,i);
    1370                 :            :   /* MT(i,i) unused */
    1371                 :            : 
    1372                 :         21 :   av2 = avma;
    1373                 :            :   /* N = itos(gdiv(mpfact(n), mpfact(n >> 1))) >> (n >> 1); */
    1374                 :            :   /* n = 2k = 12; N = (2k)! / (k! * 2^k) = 10395 */
    1375                 :         21 :   N = 10395;
    1376         [ -  + ]:         21 :   if (DEBUGLEVEL>=4) err_printf("A4GaloisConj:will test %ld permutations\n", N);
    1377                 :         21 :   uel(ar,4) = umael(MT,11,12);
    1378                 :         21 :   uel(ar,3) = uel(ar,4) + umael(MT,9,10);
    1379                 :         21 :   uel(ar,2) = uel(ar,3) + umael(MT,7,8);
    1380                 :         21 :   uel(ar,1) = uel(ar,2) + umael(MT,5,6);
    1381         [ +  - ]:      61159 :   for (i = 0; i < N; i++)
    1382                 :            :   {
    1383                 :            :     long g;
    1384         [ +  + ]:      61159 :     if (i)
    1385                 :            :     {
    1386                 :      61138 :       long a, x = i, y = 1;
    1387         [ +  + ]:      86198 :       do { y += 2; a = x%y; x = x/y; } while (!a);
    1388   [ +  +  +  +  :      61138 :       switch (y)
                   +  - ]
    1389                 :            :       {
    1390                 :            :       case 3:
    1391                 :      40768 :         lswap(t[2], t[2-a]);
    1392                 :      40768 :         break;
    1393                 :            :       case 5:
    1394                 :      16296 :         x = t[0]; t[0] = t[2]; t[2] = t[1]; t[1] = x;
    1395                 :      16296 :         lswap(t[4], t[4-a]);
    1396                 :      16296 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1397                 :      16296 :         break;
    1398                 :            :       case 7:
    1399                 :       3507 :         x = t[0]; t[0] = t[4]; t[4] = t[3]; t[3] = t[1]; t[1] = t[2]; t[2] = x;
    1400                 :       3507 :         lswap(t[6], t[6-a]);
    1401                 :       3507 :         uel(ar,2) = uel(ar,3) + umael(MT,t[6],t[7]);
    1402                 :       3507 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1403                 :       3507 :         break;
    1404                 :            :       case 9:
    1405                 :        518 :         x = t[0]; t[0] = t[6]; t[6] = t[5]; t[5] = t[3]; t[3] = x;
    1406                 :        518 :         lswap(t[1], t[4]);
    1407                 :        518 :         lswap(t[8], t[8-a]);
    1408                 :        518 :         uel(ar,3) = uel(ar,4) + umael(MT,t[8],t[9]);
    1409                 :        518 :         uel(ar,2) = uel(ar,3) + umael(MT,t[6],t[7]);
    1410                 :        518 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1411                 :        518 :         break;
    1412                 :            :       case 11:
    1413                 :         49 :         x = t[0]; t[0] = t[8]; t[8] = t[7]; t[7] = t[5]; t[5] = t[1];
    1414                 :         49 :         t[1] = t[6]; t[6] = t[3]; t[3] = t[2]; t[2] = t[4]; t[4] = x;
    1415                 :         49 :         lswap(t[10], t[10-a]);
    1416                 :         49 :         uel(ar,4) = umael(MT,t[10],t[11]);
    1417                 :         49 :         uel(ar,3) = uel(ar,4) + umael(MT,t[8],t[9]);
    1418                 :         49 :         uel(ar,2) = uel(ar,3) + umael(MT,t[6],t[7]);
    1419                 :         49 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1420                 :            :       }
    1421                 :            :     }
    1422                 :      61159 :     g = uel(ar,1)+umael(MT,t[0],t[1])+umael(MT,t[2],t[3]);
    1423         [ +  + ]:      61159 :     if (headlongisint(g,n))
    1424                 :            :     {
    1425         [ +  + ]:        147 :       for (k = 0; k < n; k += 2)
    1426                 :            :       {
    1427                 :        126 :         pft[t[k]] = t[k+1];
    1428                 :        126 :         pft[t[k+1]] = t[k];
    1429                 :            :       }
    1430         [ +  - ]:         21 :       if (galois_test_perm(td, pft)) break;
    1431                 :          0 :       hop++;
    1432                 :            :     }
    1433                 :      61138 :     avma = av2;
    1434                 :            :   }
    1435 [ -  + ][ #  # ]:         21 :   if (DEBUGLEVEL >= 1 && hop)
    1436                 :          0 :     err_printf("A4GaloisConj: %ld hop over %ld iterations\n", hop, N);
    1437         [ -  + ]:         21 :   if (i == N) { avma = ltop; return gen_0; }
    1438                 :            :   /* N = itos(gdiv(mpfact(n >> 1), mpfact(n >> 2))) >> 1; */
    1439                 :         21 :   N = 60;
    1440         [ -  + ]:         21 :   if (DEBUGLEVEL >= 4) err_printf("A4GaloisConj:sigma=%Ps \n", pft);
    1441         [ +  + ]:         84 :   for (k = 0; k < n; k += 4)
    1442                 :            :   {
    1443                 :         63 :     u[k+3] = t[k+3];
    1444                 :         63 :     u[k+2] = t[k+1];
    1445                 :         63 :     u[k+1] = t[k+2];
    1446                 :         63 :     u[k]   = t[k];
    1447                 :            :   }
    1448         [ +  - ]:        973 :   for (i = 0; i < N; i++)
    1449                 :            :   {
    1450                 :        973 :     ulong g = 0;
    1451         [ +  + ]:        973 :     if (i)
    1452                 :            :     {
    1453                 :        952 :       long a, x = i, y = -2;
    1454         [ +  + ]:       1505 :       do { y += 4; a = x%y; x = x/y; } while (!a);
    1455                 :        952 :       lswap(u[0],u[2]);
    1456   [ +  +  +  - ]:        952 :       switch (y)
    1457                 :            :       {
    1458                 :            :       case 2:
    1459                 :        476 :         break;
    1460                 :            :       case 6:
    1461                 :        399 :         lswap(u[4],u[6]);
    1462         [ +  + ]:        399 :         if (!(a & 1))
    1463                 :            :         {
    1464                 :        161 :           a = 4 - (a>>1);
    1465                 :        161 :           lswap(u[6], u[a]);
    1466                 :        161 :           lswap(u[4], u[a-2]);
    1467                 :            :         }
    1468                 :        399 :         break;
    1469                 :            :       case 10:
    1470                 :         77 :         x = u[6];
    1471                 :         77 :         u[6] = u[3];
    1472                 :         77 :         u[3] = u[2];
    1473                 :         77 :         u[2] = u[4];
    1474                 :         77 :         u[4] = u[1];
    1475                 :         77 :         u[1] = u[0];
    1476                 :         77 :         u[0] = x;
    1477         [ +  + ]:         77 :         if (a >= 3) a += 2;
    1478                 :         77 :         a = 8 - a;
    1479                 :         77 :         lswap(u[10],u[a]);
    1480                 :         77 :         lswap(u[8], u[a-2]);
    1481                 :         77 :         break;
    1482                 :            :       }
    1483                 :            :     }
    1484         [ +  + ]:       6811 :     for (k = 0; k < n; k += 2) g += mael(MT,u[k],u[k+1]);
    1485         [ +  + ]:        973 :     if (headlongisint(g,n))
    1486                 :            :     {
    1487         [ +  + ]:        147 :       for (k = 0; k < n; k += 2)
    1488                 :            :       {
    1489                 :        126 :         pfu[u[k]] = u[k+1];
    1490                 :        126 :         pfu[u[k+1]] = u[k];
    1491                 :            :       }
    1492         [ +  - ]:         21 :       if (galois_test_perm(td, pfu)) break;
    1493                 :          0 :       hop++;
    1494                 :            :     }
    1495                 :        952 :     avma = av2;
    1496                 :            :   }
    1497         [ -  + ]:         21 :   if (i == N) { avma = ltop; return gen_0; }
    1498 [ -  + ][ #  # ]:         21 :   if (DEBUGLEVEL >= 1 && hop)
    1499                 :          0 :     err_printf("A4GaloisConj: %ld hop over %ld iterations\n", hop, N);
    1500         [ -  + ]:         21 :   if (DEBUGLEVEL >= 4) err_printf("A4GaloisConj:tau=%Ps \n", pfu);
    1501                 :         21 :   avma = av2;
    1502                 :         21 :   orb = mkvec2(pft,pfu);
    1503                 :         21 :   O = vecperm_orbits(orb, 12);
    1504         [ -  + ]:         21 :   if (DEBUGLEVEL >= 4) {
    1505                 :          0 :     err_printf("A4GaloisConj:orb=%Ps\n", orb);
    1506                 :          0 :     err_printf("A4GaloisConj:O=%Ps \n", O);
    1507                 :            :   }
    1508                 :         21 :   av2 = avma;
    1509                 :         21 :   O1 = gel(O,1); O2 = gel(O,2); O3 = gel(O,3);
    1510         [ +  - ]:         35 :   for (j = 0; j < 2; j++)
    1511                 :            :   {
    1512                 :         35 :     pfv[O1[1]] = O2[1];
    1513                 :         35 :     pfv[O1[2]] = O2[3+j];
    1514                 :         35 :     pfv[O1[3]] = O2[4 - (j << 1)];
    1515                 :         35 :     pfv[O1[4]] = O2[2+j];
    1516         [ +  + ]:        119 :     for (i = 0; i < 4; i++)
    1517                 :            :     {
    1518                 :        105 :       ulong g = 0;
    1519   [ +  +  +  +  :        105 :       switch (i)
                      - ]
    1520                 :            :       {
    1521                 :         35 :       case 0: break;
    1522                 :         35 :       case 1: lswap(O3[1], O3[2]); lswap(O3[3], O3[4]); break;
    1523                 :         21 :       case 2: lswap(O3[1], O3[4]); lswap(O3[2], O3[3]); break;
    1524                 :         14 :       case 3: lswap(O3[1], O3[2]); lswap(O3[3], O3[4]); break;
    1525                 :            :       }
    1526                 :        105 :       pfv[O2[1]]          = O3[1];
    1527                 :        105 :       pfv[O2[3+j]]        = O3[4-j];
    1528                 :        105 :       pfv[O2[4 - (j<<1)]] = O3[2 + (j<<1)];
    1529                 :        105 :       pfv[O2[2+j]]        = O3[3-j];
    1530                 :        105 :       pfv[O3[1]]          = O1[1];
    1531                 :        105 :       pfv[O3[4-j]]        = O1[2];
    1532                 :        105 :       pfv[O3[2 + (j<<1)]] = O1[3];
    1533                 :        105 :       pfv[O3[3-j]]        = O1[4];
    1534         [ +  + ]:       1365 :       for (k = 1; k <= n; k++) g += mael(mt,k,pfv[k]);
    1535 [ +  + ][ +  - ]:        105 :       if (headlongisint(g,n) && galois_test_perm(td, pfv))
    1536                 :            :       {
    1537                 :         21 :         avma = av;
    1538         [ -  + ]:         21 :         if (DEBUGLEVEL >= 1)
    1539                 :          0 :           err_printf("A4GaloisConj:%ld hop over %d iterations max\n",
    1540                 :            :                      hop, 10395 + 68);
    1541                 :         21 :         return res;
    1542                 :            :       }
    1543                 :         84 :       hop++; avma = av2;
    1544                 :            :     }
    1545                 :            :   }
    1546                 :         21 :   avma = ltop; return gen_0; /* Fail */
    1547                 :            : }
    1548                 :            : 
    1549                 :            : /* S4 */
    1550                 :            : static void
    1551                 :        133 : s4makelift(GEN u, struct galois_lift *gl, GEN liftpow)
    1552                 :            : {
    1553                 :        133 :   GEN s = automorphismlift(u, gl, NULL);
    1554                 :            :   long i;
    1555                 :        133 :   gel(liftpow,1) = s;
    1556         [ +  + ]:       3059 :   for (i = 2; i < lg(liftpow); i++)
    1557                 :       2926 :     gel(liftpow,i) = FpXQ_mul(gel(liftpow,i-1), s, gl->TQ, gl->Q);
    1558                 :        133 : }
    1559                 :            : static long
    1560                 :       2338 : s4test(GEN u, GEN liftpow, struct galois_lift *gl, GEN phi)
    1561                 :            : {
    1562                 :       2338 :   pari_sp av = avma;
    1563                 :            :   GEN res, Q, Q2;
    1564                 :       2338 :   long bl, i, d = lg(u)-2;
    1565                 :            :   pari_timer ti;
    1566         [ -  + ]:       2338 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1567         [ -  + ]:       2338 :   if (!d) return 0;
    1568                 :       2338 :   Q = gl->Q; Q2 = shifti(Q,-1);
    1569                 :       2338 :   res = gel(u,2);
    1570         [ +  + ]:      54836 :   for (i = 1; i < d; i++)
    1571         [ +  - ]:      52498 :     if (lg(gel(liftpow,i))>2)
    1572                 :      52498 :       res = addii(res, mulii(gmael(liftpow,i,2), gel(u,i+2)));
    1573                 :       2338 :   res = remii(res,Q);
    1574         [ +  - ]:       2338 :   if (gl->den != gen_1) res = mulii(res, gl->den);
    1575                 :       2338 :   res = centermodii(res, Q,Q2);
    1576         [ +  + ]:       2338 :   if (absi_cmp(res, gl->gb->bornesol) > 0) { avma = av; return 0; }
    1577                 :        126 :   res = scalar_ZX_shallow(gel(u,2),varn(u));
    1578         [ +  + ]:       2562 :   for (i = 1; i < d ; i++)
    1579         [ +  - ]:       2436 :     if (lg(gel(liftpow,i))>2)
    1580                 :       2436 :       res = ZX_add(res, ZX_Z_mul(gel(liftpow,i), gel(u,i+2)));
    1581                 :        126 :   res = FpX_red(res, Q);
    1582         [ +  - ]:        126 :   if (gl->den != gen_1) res = FpX_Fp_mul(res, gl->den, Q);
    1583                 :        126 :   res = FpX_center(res, Q, shifti(Q,-1));
    1584                 :        126 :   bl = poltopermtest(res, gl, phi);
    1585         [ -  + ]:        126 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "s4test()");
    1586                 :       2338 :   avma = av; return bl;
    1587                 :            : }
    1588                 :            : 
    1589                 :            : static GEN
    1590                 :        399 : aux(long a, long b, GEN T, GEN M, GEN p, GEN *pu)
    1591                 :            : {
    1592                 :        399 :   *pu = FpX_mul(gel(T,b), gel(T,a),p);
    1593                 :        399 :   return FpX_chinese_coprime(gmael(M,a,b), gmael(M,b,a),
    1594                 :        798 :                              gel(T,b), gel(T,a), *pu, p);
    1595                 :            : }
    1596                 :            : 
    1597                 :            : static GEN
    1598                 :        133 : s4releveauto(GEN misom,GEN Tmod,GEN Tp,GEN p,long a1,long a2,long a3,long a4,long a5,long a6)
    1599                 :            : {
    1600                 :        133 :   pari_sp av = avma;
    1601                 :            :   GEN u4,u5;
    1602                 :            :   GEN pu1, pu2, pu3, pu4;
    1603                 :        133 :   GEN u1 = aux(a1, a2, Tmod, misom, p, &pu1);
    1604                 :        133 :   GEN u2 = aux(a3, a4, Tmod, misom, p, &pu2);
    1605                 :        133 :   GEN u3 = aux(a5, a6, Tmod, misom, p, &pu3);
    1606                 :        133 :   pu4 = FpX_mul(pu1,pu2,p);
    1607                 :        133 :   u4 = FpX_chinese_coprime(u1,u2,pu1,pu2,pu4,p);
    1608                 :        133 :   u5 = FpX_chinese_coprime(u4,u3,pu4,pu3,Tp,p);
    1609                 :        133 :   return gerepileupto(av, u5);
    1610                 :            : }
    1611                 :            : static GEN
    1612                 :       3773 : lincomb(GEN A, GEN B, GEN pauto, long j)
    1613                 :            : {
    1614                 :       3773 :   long k = (-j) & 3;
    1615         [ +  + ]:       3773 :   if (j == k) return ZX_mul(ZX_add(A,B), gel(pauto, j+1));
    1616                 :       3773 :   return ZX_add(ZX_mul(A, gel(pauto, j+1)), ZX_mul(B, gel(pauto, k+1)));
    1617                 :            : }
    1618                 :            : /* FIXME: could use the intheadlong technique */
    1619                 :            : static GEN
    1620                 :         21 : s4galoisgen(struct galois_lift *gl)
    1621                 :            : {
    1622                 :         21 :   const long n = 24;
    1623                 :            :   struct galois_testlift gt;
    1624                 :         21 :   pari_sp av, ltop2, ltop = avma;
    1625                 :            :   long i, j;
    1626                 :         21 :   GEN sigma, tau, phi, res, r1,r2,r3,r4, pj, p = gl->p, Q = gl->Q, TQ = gl->TQ;
    1627                 :            :   GEN sg, Tp, Tmod, isom, isominv, misom, Bcoeff, pauto, liftpow, aut;
    1628                 :            : 
    1629                 :         21 :   res = cgetg(3, t_VEC);
    1630                 :         21 :   r1 = cgetg(n+1, t_VECSMALL);
    1631                 :         21 :   r2 = cgetg(n+1, t_VECSMALL);
    1632                 :         21 :   r3 = cgetg(n+1, t_VECSMALL);
    1633                 :         21 :   r4 = cgetg(n+1, t_VECSMALL);
    1634                 :         21 :   gel(res,1)= mkvec4(r1,r2,r3,r4);
    1635                 :         21 :   gel(res,2) = mkvecsmall4(2,2,3,2);
    1636                 :         21 :   ltop2 = avma;
    1637                 :         21 :   sg = identity_perm(6);
    1638                 :         21 :   pj = zero_zv(6);
    1639                 :         21 :   sigma = cgetg(n+1, t_VECSMALL);
    1640                 :         21 :   tau = cgetg(n+1, t_VECSMALL);
    1641                 :         21 :   phi = cgetg(n+1, t_VECSMALL);
    1642                 :         21 :   Tp = FpX_red(gl->T,p);
    1643                 :         21 :   Tmod = gel(FpX_factor(Tp,p), 1);
    1644                 :         21 :   isom    = cgetg(lg(Tmod), t_VEC);
    1645                 :         21 :   isominv = cgetg(lg(Tmod), t_VEC);
    1646                 :         21 :   misom   = cgetg(lg(Tmod), t_MAT);
    1647                 :         21 :   aut = galoisdolift(gl, NULL);
    1648                 :         21 :   inittestlift(aut, Tmod, gl, &gt);
    1649                 :         21 :   Bcoeff = gt.bezoutcoeff;
    1650                 :         21 :   pauto = gt.pauto;
    1651         [ +  + ]:        147 :   for (i = 1; i < lg(isom); i++)
    1652                 :            :   {
    1653                 :        126 :     gel(misom,i) = cgetg(lg(Tmod), t_COL);
    1654                 :        126 :     gel(isom,i) = FpX_ffisom(gel(Tmod,1), gel(Tmod,i), p);
    1655         [ -  + ]:        126 :     if (DEBUGLEVEL >= 6)
    1656                 :          0 :       err_printf("S4GaloisConj:Computing isomorphisms %d:%Ps\n", i,
    1657                 :          0 :                  gel(isom,i));
    1658                 :        126 :     gel(isominv,i) = FpXQ_ffisom_inv(gel(isom,i), gel(Tmod,i),p);
    1659                 :            :   }
    1660         [ +  + ]:        147 :   for (i = 1; i < lg(isom); i++)
    1661         [ +  + ]:        882 :     for (j = 1; j < lg(isom); j++)
    1662                 :        756 :       gmael(misom,i,j) = FpX_FpXQ_eval(gel(isominv,i),gel(isom,j),
    1663                 :        756 :                                          gel(Tmod,j),p);
    1664                 :         21 :   liftpow = cgetg(24, t_VEC);
    1665                 :         21 :   av = avma;
    1666         [ +  - ]:         35 :   for (i = 0; i < 3; i++, avma = av)
    1667                 :            :   {
    1668                 :            :     pari_sp av1, av2, av3;
    1669                 :            :     GEN u, u1, u2, u3;
    1670                 :            :     long j1, j2, j3;
    1671         [ +  + ]:         35 :     if (i)
    1672                 :            :     {
    1673         [ +  - ]:         14 :       if (i == 1) { lswap(sg[2],sg[3]); }
    1674                 :          0 :       else        { lswap(sg[1],sg[3]); }
    1675                 :            :     }
    1676                 :         35 :     u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[2],sg[3],sg[4],sg[5],sg[6]);
    1677                 :         35 :     s4makelift(u, gl, liftpow);
    1678                 :         35 :     av1 = avma;
    1679         [ +  + ]:        133 :     for (j1 = 0; j1 < 4; j1++, avma = av1)
    1680                 :            :     {
    1681                 :        119 :       u1 = lincomb(gel(Bcoeff,sg[5]),gel(Bcoeff,sg[6]), pauto,j1);
    1682                 :        119 :       u1 = FpX_rem(u1, TQ, Q); av2 = avma;
    1683         [ +  + ]:        526 :       for (j2 = 0; j2 < 4; j2++, avma = av2)
    1684                 :            :       {
    1685                 :        428 :         u2 = lincomb(gel(Bcoeff,sg[3]),gel(Bcoeff,sg[4]), pauto,j2);
    1686                 :        428 :         u2 = FpX_rem(FpX_add(u1, u2, Q), TQ,Q); av3 = avma;
    1687         [ +  + ]:       2093 :         for (j3 = 0; j3 < 4; j3++, avma = av3)
    1688                 :            :         {
    1689                 :       1686 :           u3 = lincomb(gel(Bcoeff,sg[1]),gel(Bcoeff,sg[2]), pauto,j3);
    1690                 :       1686 :           u3 = FpX_rem(FpX_add(u2, u3, Q), TQ,Q);
    1691         [ -  + ]:       1686 :           if (DEBUGLEVEL >= 4)
    1692                 :          0 :             err_printf("S4GaloisConj:Testing %d/3:%d/4:%d/4:%d/4:%Ps\n",
    1693                 :            :                        i,j1,j2,j3, sg);
    1694         [ +  + ]:       1686 :           if (s4test(u3, liftpow, gl, sigma))
    1695                 :            :           {
    1696                 :         21 :             pj[1] = j3;
    1697                 :         21 :             pj[2] = j2;
    1698                 :         21 :             pj[3] = j1; goto suites4;
    1699                 :            :           }
    1700                 :            :         }
    1701                 :            :       }
    1702                 :            :     }
    1703                 :            :   }
    1704                 :          0 :   avma = ltop; return gen_0;
    1705                 :            : suites4:
    1706         [ -  + ]:         21 :   if (DEBUGLEVEL >= 4) err_printf("S4GaloisConj:sigma=%Ps\n", sigma);
    1707         [ -  + ]:         21 :   if (DEBUGLEVEL >= 4) err_printf("S4GaloisConj:pj=%Ps\n", pj);
    1708                 :         21 :   avma = av;
    1709         [ +  - ]:         42 :   for (j = 1; j <= 3; j++)
    1710                 :            :   {
    1711                 :            :     pari_sp av2, av3;
    1712                 :            :     GEN u;
    1713                 :            :     long w, l, z;
    1714                 :         42 :     z = sg[1]; sg[1] = sg[3]; sg[3] = sg[5]; sg[5] = z;
    1715                 :         42 :     z = sg[2]; sg[2] = sg[4]; sg[4] = sg[6]; sg[6] = z;
    1716                 :         42 :     z = pj[1]; pj[1] = pj[2]; pj[2] = pj[3]; pj[3] = z;
    1717         [ +  + ]:         98 :     for (l = 0; l < 2; l++, avma = av)
    1718                 :            :     {
    1719                 :         77 :       u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[3],sg[2],sg[4],sg[5],sg[6]);
    1720                 :         77 :       s4makelift(u, gl, liftpow);
    1721                 :         77 :       av2 = avma;
    1722         [ +  + ]:        210 :       for (w = 0; w < 4; w += 2, avma = av2)
    1723                 :            :       {
    1724                 :            :         GEN uu;
    1725                 :        154 :         pj[6] = (w + pj[3]) & 3;
    1726                 :        154 :         uu = lincomb(gel(Bcoeff,sg[5]),gel(Bcoeff,sg[6]), pauto, pj[6]);
    1727                 :        154 :         uu = FpX_rem(FpX_red(uu,Q), TQ, Q);
    1728                 :        154 :         av3 = avma;
    1729         [ +  + ]:        703 :         for (i = 0; i < 4; i++, avma = av3)
    1730                 :            :         {
    1731                 :            :           GEN u;
    1732                 :        570 :           pj[4] = i;
    1733                 :        570 :           pj[5] = (i + pj[2] - pj[1]) & 3;
    1734         [ -  + ]:        570 :           if (DEBUGLEVEL >= 4)
    1735                 :          0 :             err_printf("S4GaloisConj:Testing %d/3:%d/2:%d/2:%d/4:%Ps:%Ps\n",
    1736                 :            :                        j-1, w >> 1, l, i, sg, pj);
    1737                 :        570 :           u = ZX_add(lincomb(gel(Bcoeff,sg[1]),gel(Bcoeff,sg[3]), pauto,pj[4]),
    1738                 :       1710 :                      lincomb(gel(Bcoeff,sg[2]),gel(Bcoeff,sg[4]), pauto,pj[5]));
    1739                 :        570 :           u = FpX_rem(FpX_add(uu,u,Q), TQ, Q);
    1740         [ +  + ]:        570 :           if (s4test(u, liftpow, gl, tau)) goto suites4_2;
    1741                 :            :         }
    1742                 :            :       }
    1743                 :         56 :       lswap(sg[3], sg[4]);
    1744                 :         56 :       pj[2] = (-pj[2]) & 3;
    1745                 :            :     }
    1746                 :            :   }
    1747                 :          0 :   avma = ltop; return gen_0;
    1748                 :            : suites4_2:
    1749                 :         21 :   avma = av;
    1750                 :            :   {
    1751                 :         21 :     long abc = (pj[1] + pj[2] + pj[3]) & 3;
    1752                 :         21 :     long abcdef = ((abc + pj[4] + pj[5] - pj[6]) & 3) >> 1;
    1753                 :            :     GEN u;
    1754                 :            :     pari_sp av2;
    1755                 :         21 :     u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[4],sg[2],sg[5],sg[3],sg[6]);
    1756                 :         21 :     s4makelift(u, gl, liftpow);
    1757                 :         21 :     av2 = avma;
    1758         [ +  - ]:         82 :     for (j = 0; j < 8; j++)
    1759                 :            :     {
    1760                 :            :       long h, g, i;
    1761                 :         82 :       h = j & 3;
    1762                 :         82 :       g = (abcdef + ((j & 4) >> 1)) & 3;
    1763                 :         82 :       i = (h + abc - g) & 3;
    1764                 :         82 :       u = ZX_add(   lincomb(gel(Bcoeff,sg[1]), gel(Bcoeff,sg[4]), pauto, g),
    1765                 :        164 :                     lincomb(gel(Bcoeff,sg[2]), gel(Bcoeff,sg[5]), pauto, h));
    1766                 :         82 :       u = FpX_add(u, lincomb(gel(Bcoeff,sg[3]), gel(Bcoeff,sg[6]), pauto, i),Q);
    1767                 :         82 :       u = FpX_rem(u, TQ, Q);
    1768         [ -  + ]:         82 :       if (DEBUGLEVEL >= 4)
    1769                 :          0 :         err_printf("S4GaloisConj:Testing %d/8 %d:%d:%d\n", j,g,h,i);
    1770         [ +  + ]:         82 :       if (s4test(u, liftpow, gl, phi)) break;
    1771                 :         61 :       avma = av2;
    1772                 :            :     }
    1773                 :            :   }
    1774         [ -  + ]:         21 :   if (j == 8) { avma = ltop; return gen_0; }
    1775         [ +  + ]:        525 :   for (i = 1; i <= n; i++)
    1776                 :            :   {
    1777                 :        504 :     r1[i] = sigma[tau[i]];
    1778                 :        504 :     r2[i] = phi[sigma[tau[phi[i]]]];
    1779                 :        504 :     r3[i] = phi[sigma[i]];
    1780                 :        504 :     r4[i] = sigma[i];
    1781                 :            :   }
    1782                 :         21 :   avma = ltop2; return res;
    1783                 :            : }
    1784                 :            : 
    1785                 :            : static GEN
    1786                 :        329 : galoisfindgroups(GEN lo, GEN sg, long f)
    1787                 :            : {
    1788                 :        329 :   pari_sp ltop = avma;
    1789                 :            :   long i, j, k;
    1790                 :        329 :   GEN V = cgetg(lg(lo), t_VEC);
    1791         [ +  + ]:       1274 :   for(j=1,i=1; i<lg(lo); i++)
    1792                 :            :   {
    1793                 :        945 :     pari_sp av = avma;
    1794                 :        945 :     GEN loi = gel(lo,i), W = cgetg(lg(loi),t_VECSMALL);
    1795         [ +  + ]:       2681 :     for (k=1; k<lg(loi); k++) W[k] = loi[k] % f;
    1796                 :        945 :     W = vecsmall_uniq(W);
    1797         [ +  + ]:        945 :     if (zv_equal(W, sg)) gel(V,j++) = loi;
    1798                 :        945 :     avma = av;
    1799                 :            :   }
    1800                 :        329 :   setlg(V,j); return gerepilecopy(ltop,V);
    1801                 :            : }
    1802                 :            : 
    1803                 :            : static long
    1804                 :       1613 : galoisfrobeniustest(GEN aut, struct galois_lift *gl, GEN frob)
    1805                 :            : {
    1806                 :       1613 :   pari_sp av = avma;
    1807                 :       1613 :   GEN tlift = aut;
    1808                 :            :   long res;
    1809         [ +  - ]:       1613 :   if (gl->den != gen_1) tlift = FpX_Fp_mul(tlift, gl->den, gl->Q);
    1810                 :       1613 :   tlift = FpX_center(tlift, gl->Q, shifti(gl->Q,-1));
    1811                 :       1613 :   res = poltopermtest(tlift, gl, frob);
    1812                 :       1613 :   avma = av; return res;
    1813                 :            : }
    1814                 :            : 
    1815                 :            : static GEN
    1816                 :        441 : galoismakepsi(long g, GEN sg, GEN pf)
    1817                 :            : {
    1818                 :        441 :   GEN psi=cgetg(g+1,t_VECSMALL);
    1819                 :            :   long i;
    1820         [ +  + ]:       1456 :   for (i = 1; i < g; i++) psi[i] = sg[pf[i]];
    1821                 :        441 :   psi[g] = sg[1]; return psi;
    1822                 :            : }
    1823                 :            : 
    1824                 :            : static GEN
    1825                 :       1624 : galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden,
    1826                 :            :     struct galois_frobenius *gf,  struct galois_borne *gb)
    1827                 :            : {
    1828                 :       1624 :   pari_sp ltop=avma, av2;
    1829                 :            :   struct galois_testlift gt;
    1830                 :            :   struct galois_lift gl;
    1831                 :       1624 :   long i, j, k, n = lg(L)-1, deg = 1, g = lg(gf->Tmod)-1;
    1832                 :       1624 :   GEN F,Fp,Fe, aut, frob, ip = utoipos(gf->p), res = cgetg(lg(L), t_VECSMALL);
    1833         [ -  + ]:       1624 :   if (DEBUGLEVEL >= 4)
    1834                 :          0 :     err_printf("GaloisConj:p=%ld deg=%ld fp=%ld\n", gf->p, deg, gf->fp);
    1835                 :       1624 :   gf->psi = const_vecsmall(g,1);
    1836                 :       1624 :   av2 = avma;
    1837                 :       1624 :   initlift(T, den, ip, L, Lden, gb, &gl);
    1838                 :       1624 :   aut = galoisdolift(&gl, res);
    1839 [ +  + ][ +  + ]:       1624 :   if (!aut || galoisfrobeniustest(aut,&gl,res))
    1840                 :            :   {
    1841                 :       1106 :     avma = av2; gf->deg = gf->fp; return res;
    1842                 :            :   }
    1843                 :        518 :   inittestlift(aut,gf->Tmod, &gl, &gt);
    1844                 :        518 :   gt.C = cgetg(gf->fp+1,t_VEC);
    1845                 :        518 :   gt.Cd= cgetg(gf->fp+1,t_VEC);
    1846         [ +  + ]:       3906 :   for (i = 1; i <= gf->fp; i++) {
    1847                 :       3388 :     gel(gt.C,i)  = zero_zv(gt.g);
    1848                 :       3388 :     gel(gt.Cd,i) = zero_zv(gt.g);
    1849                 :            :   }
    1850                 :            : 
    1851                 :        518 :   F =factoru(gf->fp);
    1852                 :        518 :   Fp = gel(F,1);
    1853                 :        518 :   Fe = gel(F,2);
    1854                 :        518 :   frob = cgetg(lg(L), t_VECSMALL);
    1855         [ +  + ]:       1127 :   for(k=lg(Fp)-1;k>=1;k--)
    1856                 :            :   {
    1857                 :        609 :     pari_sp btop=avma;
    1858                 :        609 :     GEN psi=NULL, fres=NULL, sg = identity_perm(1);
    1859                 :        609 :     long el=gf->fp, dg=1, dgf=1, e, pr;
    1860         [ +  + ]:       1386 :     for(e=1; e<=Fe[k]; e++)
    1861                 :            :     {
    1862                 :            :       GEN lo, pf;
    1863                 :            :       long l;
    1864                 :        994 :       dg *= Fp[k]; el /= Fp[k];
    1865         [ -  + ]:        994 :       if (DEBUGLEVEL>=4) err_printf("Trying degre %d.\n",dg);
    1866         [ +  + ]:        994 :       if (galoisfrobeniustest(gel(gt.pauto,el+1),&gl,frob))
    1867                 :            :       {
    1868                 :        336 :         psi = const_vecsmall(g,1);
    1869                 :        336 :         dgf = dg; fres = gcopy(frob); continue;
    1870                 :            :       }
    1871                 :        658 :       lo = listznstarelts(dg, n / gf->fp);
    1872         [ +  + ]:        658 :       if (e!=1) lo = galoisfindgroups(lo, sg, dgf);
    1873         [ -  + ]:        658 :       if (DEBUGLEVEL>=4) err_printf("Galoisconj:Subgroups list:%Ps\n", lo);
    1874         [ +  + ]:       1372 :       for (l = 1; l < lg(lo); l++)
    1875 [ +  + ][ +  + ]:       1155 :         if (lg(gel(lo,l))>2 && frobeniusliftall(gel(lo,l), el, &pf,&gl,&gt, frob))
    1876                 :            :         {
    1877                 :        441 :           sg  = gcopy(gel(lo,l));
    1878                 :        441 :           psi = galoismakepsi(g,sg,pf);
    1879                 :        441 :           dgf = dg; fres = gcopy(frob); break;
    1880                 :            :         }
    1881         [ +  + ]:        658 :       if (l == lg(lo)) break;
    1882                 :            :     }
    1883         [ +  + ]:        609 :     if (dgf == 1) { avma = btop; continue; }
    1884                 :        476 :     pr = deg*dgf;
    1885         [ +  + ]:        476 :     if (deg == 1)
    1886                 :            :     {
    1887         [ +  + ]:      11620 :       for(i=1;i<lg(res);i++) res[i]=fres[i];
    1888         [ +  + ]:       2303 :       for(i=1;i<lg(psi);i++) gf->psi[i]=psi[i];
    1889                 :            :     }
    1890                 :            :     else
    1891                 :            :     {
    1892                 :         56 :       GEN cp = perm_mul(res,fres);
    1893         [ +  + ]:       1694 :       for(i=1;i<lg(res);i++) res[i] = cp[i];
    1894         [ +  + ]:        210 :       for(i=1;i<lg(psi);i++) gf->psi[i] = (dgf*gf->psi[i] + deg*psi[i]) % pr;
    1895                 :            :     }
    1896                 :        476 :     deg = pr; avma = btop;
    1897                 :            :   }
    1898         [ +  + ]:       3906 :   for (i = 1; i <= gf->fp; i++)
    1899         [ +  + ]:      17969 :     for (j = 1; j <= gt.g; j++)
    1900         [ +  + ]:      14581 :       if (mael(gt.C,i,j)) gunclone(gmael(gt.C,i,j));
    1901 [ -  + ][ #  # ]:        518 :   if (DEBUGLEVEL>=4 && res) err_printf("Best lift: %d\n",deg);
    1902         [ +  + ]:        518 :   if (deg==1) { avma = ltop; return NULL; }
    1903                 :            :   else
    1904                 :            :   {
    1905                 :            :     /* Normalize result so that psi[g]=1 */
    1906                 :        420 :     long im = Fl_inv(gf->psi[g], deg);
    1907                 :        420 :     GEN cp = perm_pow(res, im);
    1908         [ +  + ]:      11620 :     for(i=1;i<lg(res);i++) res[i] = cp[i];
    1909         [ +  + ]:       2303 :     for(i=1;i<lg(gf->psi);i++) gf->psi[i] = Fl_mul(im, gf->psi[i], deg);
    1910                 :       1624 :     avma = av2; gf->deg = deg; return res;
    1911                 :            :   }
    1912                 :            : }
    1913                 :            : 
    1914                 :            : /* return NULL if not Galois */
    1915                 :            : static GEN
    1916                 :       1526 : galoisfindfrobenius(GEN T, GEN L, GEN den, struct galois_frobenius *gf,
    1917                 :            :     struct galois_borne *gb, const struct galois_analysis *ga)
    1918                 :            : {
    1919                 :       1526 :   pari_sp ltop = avma, av;
    1920         [ +  + ]:       1526 :   long Try = 0, n = degpol(T), deg, gmask = (ga->group&ga_ext_2)? 3: 1;
    1921                 :       1526 :   GEN frob, Lden = makeLden(L,den,gb);
    1922                 :            :   forprime_t S;
    1923                 :            : 
    1924                 :       1526 :   u_forprime_init(&S, ga->p, ULONG_MAX);
    1925                 :       1526 :   av = avma;
    1926                 :       1526 :   deg = gf->deg = ga->deg;
    1927         [ +  - ]:       1631 :   while ((gf->p = u_forprime_next(&S)))
    1928                 :            :   {
    1929                 :            :     pari_sp lbot;
    1930                 :            :     GEN Ti, Tp;
    1931                 :            :     long nb, d;
    1932                 :       1631 :     avma = av;
    1933                 :       1631 :     Tp = ZX_to_Flx(T, gf->p);
    1934         [ -  + ]:       1631 :     if (!Flx_is_squarefree(Tp, gf->p)) continue;
    1935                 :       1631 :     Ti = gel(Flx_factor(Tp, gf->p), 1);
    1936                 :       1631 :     nb = lg(Ti)-1; d = degpol(gel(Ti,1));
    1937 [ +  + ][ -  + ]:       1631 :     if (nb > 1 && degpol(gel(Ti,nb)) != d) { avma = ltop; return NULL; }
    1938 [ +  + ][ -  + ]:       1631 :     if (((gmask&1)==0 || d % deg) && ((gmask&2)==0 || odd(d))) continue;
         [ +  - ][ +  + ]
    1939         [ -  + ]:       1624 :     if (DEBUGLEVEL >= 1) err_printf("GaloisConj: Trying p=%ld\n", gf->p);
    1940                 :       1624 :     FlxV_to_ZXV_inplace(Ti);
    1941                 :       1624 :     gf->fp = d;
    1942                 :       1624 :     gf->Tmod = Ti; lbot = avma;
    1943                 :       1624 :     frob = galoisfrobeniuslift(T, den, L, Lden, gf, gb);
    1944         [ +  + ]:       1624 :     if (frob)
    1945                 :            :     {
    1946                 :            :       GEN *gptr[3];
    1947                 :       1526 :       gf->Tmod = gcopy(Ti);
    1948                 :       1526 :       gptr[0]=&gf->Tmod; gptr[1]=&gf->psi; gptr[2]=&frob;
    1949                 :       1526 :       gerepilemanysp(ltop,lbot,gptr,3); return frob;
    1950                 :            :     }
    1951 [ +  + ][ +  + ]:         98 :     if ((ga->group&ga_all_normal) && d % deg == 0) gmask &= ~1;
    1952                 :            :     /* The first prime degree is always divisible by deg, so we don't
    1953                 :            :      * have to worry about ext_2 being used before regular supersolvable*/
    1954         [ -  + ]:         98 :     if (!gmask) { avma = ltop; return NULL; }
    1955 [ +  + ][ -  + ]:         98 :     if ((ga->group&ga_non_wss) && ++Try > ((3*n)>>1))
    1956                 :            :     {
    1957                 :          0 :       pari_warn(warner,"Galois group probably not weakly super solvable");
    1958                 :          0 :       return NULL;
    1959                 :            :     }
    1960                 :            :   }
    1961         [ #  # ]:          0 :   if (!gf->p) pari_err_OVERFLOW("galoisfindfrobenius [ran out of primes]");
    1962                 :       1526 :   return NULL;
    1963                 :            : }
    1964                 :            : 
    1965                 :            : /* compute g such that tau(Pmod[#])= tau(Pmod[g]) */
    1966                 :            : 
    1967                 :            : static long
    1968                 :       1001 : get_image(GEN tau, GEN P, GEN Pmod, GEN p)
    1969                 :            : {
    1970                 :       1001 :   pari_sp av = avma;
    1971                 :       1001 :   long g, gp = lg(Pmod)-1;
    1972                 :       1001 :   tau = RgX_to_FpX(tau, p);
    1973                 :       1001 :   tau = FpX_FpXQ_eval(gel(Pmod, gp), tau, P, p);
    1974                 :       1001 :   tau = FpX_normalize(FpX_gcd(P, tau, p), p);
    1975         [ +  - ]:       2296 :   for (g = 1; g <= gp; g++)
    1976         [ +  + ]:       2296 :     if (ZX_equal(tau, gel(Pmod,g))) { avma = av; return g; }
    1977                 :       1001 :   avma = av; return 0;
    1978                 :            : }
    1979                 :            : 
    1980                 :            : static GEN
    1981                 :            : galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne *gb,
    1982                 :            :           const struct galois_analysis *ga);
    1983                 :            : static GEN
    1984                 :        770 : galoisgenfixedfield(GEN Tp, GEN Pmod, GEN V, GEN ip, struct galois_borne *gb)
    1985                 :            : {
    1986                 :            :   GEN  P, PL, Pden, PM, Pp;
    1987                 :            :   GEN  tau, PG, Pg;
    1988                 :            :   long g, lP;
    1989                 :        770 :   long x=varn(Tp);
    1990                 :        770 :   P=gel(V,3);
    1991                 :        770 :   PL=gel(V,2);
    1992                 :        770 :   Pp = FpX_red(P,ip);
    1993         [ -  + ]:        770 :   if (DEBUGLEVEL>=6)
    1994                 :          0 :     err_printf("GaloisConj: Fixed field %Ps\n",P);
    1995         [ +  + ]:        770 :   if (degpol(P)==2)
    1996                 :            :   {
    1997                 :        490 :     PG=cgetg(3,t_VEC);
    1998                 :        490 :     gel(PG,1) = mkvec( mkvecsmall2(2,1) );
    1999                 :        490 :     gel(PG,2) = mkvecsmall(2);
    2000                 :        490 :     tau = deg1pol_shallow(gen_m1, negi(gel(P,3)), x);
    2001                 :        490 :     g = get_image(tau, Pp, Pmod, ip);
    2002         [ -  + ]:        490 :     if (!g) return NULL;
    2003                 :        490 :     Pg = mkvecsmall(g);
    2004                 :            :   }
    2005                 :            :   else
    2006                 :            :   {
    2007                 :            :     struct galois_analysis Pga;
    2008                 :            :     struct galois_borne Pgb;
    2009                 :            :     GEN mod, mod2;
    2010                 :            :     long j;
    2011         [ +  + ]:        280 :     if (!galoisanalysis(P, &Pga, 0)) return NULL;
    2012                 :        266 :     Pgb.l = gb->l;
    2013                 :        266 :     Pden = galoisborne(P, NULL, &Pgb, degpol(P));
    2014                 :            : 
    2015         [ +  + ]:        266 :     if (Pgb.valabs > gb->valabs)
    2016                 :            :     {
    2017         [ -  + ]:         35 :       if (DEBUGLEVEL>=4)
    2018                 :          0 :         err_printf("GaloisConj:increase prec of p-adic roots of %ld.\n"
    2019                 :          0 :             ,Pgb.valabs-gb->valabs);
    2020                 :         35 :       PL = ZpX_liftroots(P,PL,gb->l,Pgb.valabs);
    2021                 :            :     }
    2022         [ +  - ]:        231 :     else if (Pgb.valabs < gb->valabs)
    2023                 :        231 :       PL = FpC_red(PL, Pgb.ladicabs);
    2024                 :        266 :     PM = vandermondeinversemod(PL, P, Pden, Pgb.ladicabs);
    2025                 :        266 :     PG = galoisgen(P, PL, PM, Pden, &Pgb, &Pga);
    2026         [ -  + ]:        266 :     if (PG == gen_0) return NULL;
    2027                 :        266 :     lP = lg(gel(PG,1));
    2028                 :        266 :     mod = Pgb.ladicabs; mod2 = shifti(mod, -1);
    2029                 :        266 :     Pg = cgetg(lP, t_VECSMALL);
    2030         [ +  + ]:        791 :     for (j = 1; j < lP; j++)
    2031                 :            :     {
    2032                 :        511 :       pari_sp btop=avma;
    2033                 :        511 :       tau = permtopol(gmael(PG,1,j), PL, PM, Pden, mod, mod2, x);
    2034                 :        511 :       g = get_image(tau, Pp, Pmod, ip);
    2035         [ -  + ]:        511 :       if (!g) return NULL;
    2036                 :        511 :       Pg[j] = g;
    2037                 :        511 :       avma = btop;
    2038                 :            :     }
    2039                 :            :   }
    2040                 :        770 :   return mkvec2(PG,Pg);
    2041                 :            : }
    2042                 :            : 
    2043                 :            : /* Let sigma^m=1,  tau*sigma*tau^-1=sigma^s.
    2044                 :            :  * Compute n so that (sigma*tau)^e = sigma^n*tau^e
    2045                 :            :  * We have n = sum_{k=0}^{e-1} s^k mod m.
    2046                 :            :  * so n*(1-s) = 1-s^e mod m
    2047                 :            :  * Unfortunately (1-s) might not invertible mod m.
    2048                 :            :  */
    2049                 :            : 
    2050                 :            : static long
    2051                 :       2408 : stpow(long s, long e, long m)
    2052                 :            : {
    2053                 :            :   long i;
    2054                 :       2408 :   long n = 1;
    2055         [ +  + ]:       3927 :   for (i = 1; i < e; i++)
    2056                 :       1519 :     n = (1 + n * s) % m;
    2057                 :       2408 :   return n;
    2058                 :            : }
    2059                 :            : 
    2060                 :            : static GEN
    2061                 :       1001 : wpow(long s, long m, long e, long n)
    2062                 :            : {
    2063                 :       1001 :   GEN   w = cgetg(n+1,t_VECSMALL);
    2064                 :       1001 :   long si = s;
    2065                 :            :   long i;
    2066                 :       1001 :   w[1] = 1;
    2067         [ +  + ]:       1204 :   for(i=2; i<=n; i++) w[i] = w[i-1]*e;
    2068         [ +  + ]:       2205 :   for(i=n; i>=1; i--)
    2069                 :            :   {
    2070                 :       1204 :     si = Fl_powu(si,e,m);
    2071                 :       1204 :     w[i] = Fl_mul(s-1, stpow(si, w[i], m), m);
    2072                 :            :   }
    2073                 :       1001 :   return w;
    2074                 :            : }
    2075                 :            : 
    2076                 :            : static GEN
    2077                 :       1001 : galoisgenliftauto(GEN O, GEN gj, long s, long n, struct galois_test *td)
    2078                 :            : {
    2079                 :       1001 :   pari_sp av = avma;
    2080                 :            :   long sr, k;
    2081                 :       1001 :   long deg = lg(gel(O,1))-1;
    2082                 :       1001 :   GEN  X  = cgetg(lg(O), t_VECSMALL);
    2083                 :       1001 :   GEN  oX = cgetg(lg(O), t_VECSMALL);
    2084                 :       1001 :   GEN  B  = perm_cycles(gj);
    2085                 :       1001 :   long oj = lg(gel(B,1)) - 1;
    2086                 :       1001 :   GEN  F  = factoru(oj);
    2087                 :       1001 :   GEN  Fp = gel(F,1);
    2088                 :       1001 :   GEN  Fe = gel(F,2);
    2089                 :       1001 :   GEN  pf = identity_perm(n);
    2090         [ -  + ]:       1001 :   if (DEBUGLEVEL >= 6)
    2091                 :          0 :     err_printf("GaloisConj: %Ps of relative order %d\n", gj, oj);
    2092         [ +  + ]:       2002 :   for (k=lg(Fp)-1; k>=1; k--)
    2093                 :            :   {
    2094                 :       1001 :     long f, dg = 1, el = oj, osel = 1, a = 0;
    2095                 :       1001 :     long p  = Fp[k], e  = Fe[k], op = oj / upowuu(p,e);
    2096                 :            :     long i;
    2097                 :       1001 :     GEN  pf1 = NULL, w, wg, Be = cgetg(e+1,t_VEC);
    2098                 :       1001 :     gel(Be,e) = cyc_pow(B, op);
    2099         [ +  + ]:       1204 :     for(i=e-1; i>=1; i--) gel(Be,i) = cyc_pow(gel(Be,i+1), p);
    2100                 :       1001 :     w = wpow(Fl_powu(s,op,deg),deg,p,e);
    2101                 :       1001 :     wg = cgetg(e+2,t_VECSMALL);
    2102                 :       1001 :     wg[e+1] = deg;
    2103         [ +  + ]:       2205 :     for (i=e; i>=1; i--) wg[i] = ugcd(wg[i+1], w[i]);
    2104         [ +  + ]:       8246 :     for (i=1; i<lg(O); i++) oX[i] = 0;
    2105         [ +  + ]:       2205 :     for (f=1; f<=e; f++)
    2106                 :            :     {
    2107                 :            :       long sel, t;
    2108                 :       1204 :       GEN Bel = gel(Be,f);
    2109                 :       1204 :       dg *= p; el /= p;
    2110                 :       1204 :       sel = Fl_powu(s,el,deg);
    2111         [ -  + ]:       1204 :       if (DEBUGLEVEL >= 6) err_printf("GaloisConj: B=%Ps\n", Bel);
    2112                 :       1204 :       sr  = cgcd(stpow(sel,p,deg),deg);
    2113         [ -  + ]:       1204 :       if (DEBUGLEVEL >= 6)
    2114                 :          0 :         err_printf("GaloisConj: exp %d: s=%ld [%ld] a=%ld w=%ld wg=%ld sr=%ld\n",
    2115                 :          0 :             dg, sel, deg, a, w[f], wg[f+1], sr);
    2116         [ +  - ]:       1365 :       for (t = 0; t < sr; t++)
    2117         [ +  + ]:       1365 :         if ((a+t*w[f])%wg[f+1]==0)
    2118                 :            :         {
    2119                 :            :           long i, j, k, st;
    2120         [ +  + ]:      13755 :           for (i = 1; i < lg(X); i++) X[i] = 0;
    2121         [ +  + ]:       6356 :           for (i = 0; i < lg(X)-1; i+=dg)
    2122         [ +  + ]:      11130 :             for (j = 1, k = p, st = t; k <= dg; j++, k += p)
    2123                 :            :             {
    2124                 :       6118 :               X[k+i] = (oX[j+i] + st)%deg;
    2125                 :       6118 :               st = (t + st*osel)%deg;
    2126                 :            :             }
    2127                 :       1344 :           pf1 = testpermutation(O, Bel, X, sel, p, sr, td);
    2128         [ +  + ]:       1344 :           if (pf1) break;
    2129                 :            :         }
    2130         [ -  + ]:       1204 :       if (!pf1) return NULL;
    2131         [ +  + ]:      11501 :       for (i=1; i<lg(O); i++) oX[i] = X[i];
    2132                 :       1204 :       osel = sel; a = (a+t*w[f])%deg;
    2133                 :            :     }
    2134                 :       1001 :     pf = perm_mul(pf, perm_pow(pf1, el));
    2135                 :            :   }
    2136                 :       1001 :   return gerepileuptoleaf(av, pf);
    2137                 :            : }
    2138                 :            : 
    2139                 :            : static GEN
    2140                 :       1568 : galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne *gb,
    2141                 :            :           const struct galois_analysis *ga)
    2142                 :            : {
    2143                 :            :   struct galois_test td;
    2144                 :            :   struct galois_frobenius gf;
    2145                 :       1568 :   pari_sp ltop = avma;
    2146                 :       1568 :   long p, deg, x, j, n = degpol(T), lP;
    2147                 :            :   GEN sigma, Tmod, res, res1, ip, frob, O, PG, PG1, PG2, Pg;
    2148                 :            : 
    2149         [ -  + ]:       1568 :   if (!ga->deg) return gen_0;
    2150                 :       1568 :   x = varn(T);
    2151         [ -  + ]:       1568 :   if (DEBUGLEVEL >= 9) err_printf("GaloisConj:denominator:%Ps\n", den);
    2152 [ +  + ][ +  + ]:       1568 :   if (n == 12 && ga->ord==3 && !ga->p4)
                 [ +  - ]
    2153                 :            :   { /* A4 is very probable: test it first */
    2154                 :         21 :     pari_sp av = avma;
    2155         [ -  + ]:         21 :     if (DEBUGLEVEL >= 4) err_printf("GaloisConj:Testing A4 first\n");
    2156                 :         21 :     inittest(L, M, gb->bornesol, gb->ladicsol, &td);
    2157                 :         21 :     PG = a4galoisgen(&td);
    2158                 :         21 :     freetest(&td);
    2159         [ +  - ]:         21 :     if (PG != gen_0) return gerepileupto(ltop, PG);
    2160                 :          0 :     avma = av;
    2161                 :            :   }
    2162 [ +  + ][ +  + ]:       1547 :   if (n == 24 && ga->ord==3)
    2163                 :            :   { /* S4 is very probable: test it first */
    2164                 :         21 :     pari_sp av = avma;
    2165                 :            :     struct galois_lift gl;
    2166         [ -  + ]:         21 :     if (DEBUGLEVEL >= 4) err_printf("GaloisConj:Testing S4 first\n");
    2167                 :         21 :     initlift(T, den, stoi(ga->p4), L, makeLden(L,den,gb), gb, &gl);
    2168                 :         21 :     PG = s4galoisgen(&gl);
    2169         [ +  - ]:         21 :     if (PG != gen_0) return gerepileupto(ltop, PG);
    2170                 :          0 :     avma = av;
    2171                 :            :   }
    2172                 :       1526 :   frob = galoisfindfrobenius(T, L, den, &gf, gb, ga);
    2173         [ -  + ]:       1526 :   if (!frob) { avma=ltop; return gen_0; }
    2174                 :       1526 :   p = gf.p; ip = utoipos(p);
    2175                 :       1526 :   Tmod = gf.Tmod;
    2176                 :       1526 :   O = perm_cycles(frob);
    2177                 :       1526 :   deg = lg(gel(O,1))-1;
    2178         [ -  + ]:       1526 :   if (DEBUGLEVEL >= 9) err_printf("GaloisConj:Orbite:%Ps\n", O);
    2179         [ +  + ]:       1526 :   if (deg == n)        /* cyclic */
    2180                 :        756 :     return gerepilecopy(ltop, mkvec2(mkvec(frob), mkvecsmall(deg)));
    2181                 :        770 :   sigma = permtopol(frob, L, M, den, gb->ladicabs, shifti(gb->ladicabs,-1), x);
    2182         [ -  + ]:        770 :   if (DEBUGLEVEL >= 9) err_printf("GaloisConj:Frobenius:%Ps\n", sigma);
    2183                 :            :   {
    2184                 :        770 :     pari_sp btop=avma;
    2185                 :            :     GEN V, Tp, Sp, Pmod;
    2186                 :        770 :     GEN OL = fixedfieldorbits(O,L);
    2187                 :        770 :     V  = fixedfieldsympol(OL, gb->ladicabs, gb->l, ip, x);
    2188                 :        770 :     Tp = FpX_red(T,ip);
    2189                 :        770 :     Sp = sympol_aut_evalmod(gel(V,1),deg,sigma,Tp,ip);
    2190                 :        770 :     Pmod = fixedfieldfactmod(Sp,ip,Tmod);
    2191                 :        770 :     PG = galoisgenfixedfield(Tp, Pmod, V, ip, gb);
    2192         [ +  + ]:        770 :     if (PG == NULL) { avma = ltop; return gen_0; }
    2193         [ -  + ]:        756 :     if (DEBUGLEVEL >= 4) err_printf("GaloisConj:Back to Earth:%Ps\n", PG);
    2194                 :        756 :     PG=gerepilecopy(btop, PG);
    2195                 :            :   }
    2196                 :        756 :   inittest(L, M, gb->bornesol, gb->ladicsol, &td);
    2197                 :        756 :   PG1 = gmael(PG, 1, 1); lP = lg(PG1);
    2198                 :        756 :   PG2 = gmael(PG, 1, 2);
    2199                 :        756 :   Pg = gel(PG, 2);
    2200                 :        756 :   res = cgetg(3, t_VEC);
    2201                 :        756 :   gel(res,1) = res1 = cgetg(lP + 1, t_VEC);
    2202                 :        756 :   gel(res,2) = vecsmall_prepend(PG2, deg);
    2203                 :        756 :   gel(res1, 1) = vecsmall_copy(frob);
    2204         [ +  + ]:       1757 :   for (j = 1; j < lP; j++)
    2205                 :            :   {
    2206                 :       1001 :     GEN pf = galoisgenliftauto(O, gel(PG1, j), gf.psi[Pg[j]], n, &td);
    2207         [ -  + ]:       1001 :     if (!pf) { freetest(&td); avma = ltop; return gen_0; }
    2208                 :       1001 :     gel(res1, j+1) = pf;
    2209                 :            :   }
    2210         [ -  + ]:        756 :   if (DEBUGLEVEL >= 4) err_printf("GaloisConj:Fini!\n");
    2211                 :        756 :   freetest(&td);
    2212                 :       1568 :   return gerepileupto(ltop, res);
    2213                 :            : }
    2214                 :            : 
    2215                 :            : /* T = polcyclo(N) */
    2216                 :            : static GEN
    2217                 :         77 : conjcyclo(GEN T, long N)
    2218                 :            : {
    2219                 :         77 :   pari_sp av = avma;
    2220                 :         77 :   long i, k = 1, d = eulerphiu(N), v = varn(T);
    2221                 :         77 :   GEN L = cgetg(d+1,t_COL);
    2222         [ +  + ]:        462 :   for (i=1; i<=N; i++)
    2223         [ +  + ]:        385 :     if (ugcd(i, N)==1)
    2224                 :            :     {
    2225                 :        196 :       GEN s = monomial(gen_1, i, v);
    2226         [ +  + ]:        196 :       if (i >= d) s = ZX_rem(s, T);
    2227                 :        196 :       gel(L,k++) = s;
    2228                 :            :     }
    2229                 :         77 :   return gerepileupto(av, gen_sort(L, (void*)&gcmp, &gen_cmp_RgX));
    2230                 :            : }
    2231                 :            : 
    2232                 :            : /* T: polynomial or nf, den multiple of common denominator of solutions or
    2233                 :            :  * NULL (unknown). If T is nf, and den unknown, use den = denom(nf.zk) */
    2234                 :            : static GEN
    2235                 :       1974 : galoisconj4_main(GEN T, GEN den, long flag)
    2236                 :            : {
    2237                 :       1974 :   pari_sp ltop = avma;
    2238                 :            :   GEN nf, G, L, M, aut;
    2239                 :            :   struct galois_analysis ga;
    2240                 :            :   struct galois_borne gb;
    2241                 :            :   long n;
    2242                 :            :   pari_timer ti;
    2243                 :            : 
    2244                 :       1974 :   T = get_nfpol(T, &nf);
    2245                 :       1974 :   n = poliscyclo(T);
    2246 [ +  + ][ +  + ]:       1974 :   if (n) return flag? galoiscyclo(n, varn(T)): conjcyclo(T, n);
    2247                 :       1855 :   n = degpol(T);
    2248         [ +  + ]:       1855 :   if (nf)
    2249         [ +  - ]:       1309 :   { if (!den) den = Q_denom(nf_get_zk(nf)); }
    2250                 :            :   else
    2251                 :            :   {
    2252         [ -  + ]:        546 :     if (n <= 0) pari_err_IRREDPOL("galoisinit",T);
    2253                 :        546 :     RgX_check_ZX(T, "galoisinit");
    2254         [ +  + ]:        546 :     if (!ZX_is_squarefree(T))
    2255                 :          7 :       pari_err_DOMAIN("galoisinit","issquarefree(pol)","=",gen_0,T);
    2256         [ +  + ]:        539 :     if (!gequal1(gel(T,n+2))) pari_err_IMPL("galoisinit(non-monic)");
    2257                 :            :   }
    2258         [ +  + ]:       1841 :   if (n == 1)
    2259                 :            :   {
    2260         [ -  + ]:          7 :     if (!flag) { G = cgetg(2, t_COL); gel(G,1) = pol_x(varn(T)); return G;}
    2261                 :          7 :     ga.l = 3;
    2262                 :          7 :     ga.deg = 1;
    2263                 :          7 :     den = gen_1;
    2264                 :            :   }
    2265         [ +  + ]:       1834 :   else if (!galoisanalysis(T, &ga, 1)) { avma = ltop; return utoipos(ga.p); }
    2266                 :            : 
    2267         [ +  + ]:       1309 :   if (den)
    2268                 :            :   {
    2269         [ -  + ]:        875 :     if (typ(den) != t_INT) pari_err_TYPE("galoisconj4 [2nd arg integer]", den);
    2270                 :        875 :     den = absi(den);
    2271                 :            :   }
    2272                 :       1309 :   gb.l = utoipos(ga.l);
    2273         [ -  + ]:       1309 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
    2274                 :       1309 :   den = galoisborne(T, den, &gb, degpol(T));
    2275         [ -  + ]:       1309 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "galoisborne()");
    2276                 :       1309 :   L = rootpadicfast(T, gb.l, gb.valabs);
    2277         [ -  + ]:       1309 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "rootpadicfast()");
    2278                 :       1309 :   M = vandermondeinversemod(L, T, den, gb.ladicabs);
    2279         [ -  + ]:       1309 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "vandermondeinversemod()");
    2280         [ +  + ]:       1309 :   if (n == 1)
    2281                 :            :   {
    2282                 :          7 :     G = cgetg(3, t_VEC);
    2283                 :          7 :     gel(G,1) = cgetg(1, t_VEC);
    2284                 :          7 :     gel(G,2) = cgetg(1, t_VECSMALL);
    2285                 :            :   }
    2286                 :            :   else
    2287                 :       1302 :     G = galoisgen(T, L, M, den, &gb, &ga);
    2288         [ -  + ]:       1309 :   if (DEBUGLEVEL >= 6) err_printf("GaloisConj:%Ps\n", G);
    2289         [ +  + ]:       1309 :   if (G == gen_0) { avma = ltop; return gen_0; }
    2290         [ -  + ]:       1295 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
    2291         [ +  + ]:       1295 :   if (flag)
    2292                 :            :   {
    2293                 :        448 :     GEN grp = cgetg(9, t_VEC);
    2294                 :        448 :     gel(grp,1) = ZX_copy(T);
    2295                 :        448 :     gel(grp,2) = mkvec3(utoipos(ga.l), utoipos(gb.valabs), icopy(gb.ladicabs));
    2296                 :        448 :     gel(grp,3) = ZC_copy(L);
    2297                 :        448 :     gel(grp,4) = ZM_copy(M);
    2298                 :        448 :     gel(grp,5) = icopy(den);
    2299                 :        448 :     gel(grp,6) = group_elts(G,n);
    2300                 :        448 :     gel(grp,7) = gcopy(gel(G,1));
    2301                 :        448 :     gel(grp,8) = gcopy(gel(G,2)); return gerepileupto(ltop, grp);
    2302                 :            :   }
    2303                 :        847 :   aut = galoisgrouptopol(group_elts(G, n),L,M,den,gb.ladicsol, varn(T));
    2304         [ -  + ]:        847 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "Computation of polynomials");
    2305                 :       1960 :   return gerepileupto(ltop, gen_sort(aut, (void*)&gcmp, &gen_cmp_RgX));
    2306                 :            : }
    2307                 :            : 
    2308                 :            : /* Heuristic computation of #Aut(T), pinit = first prime to be tested */
    2309                 :            : long
    2310                 :        455 : numberofconjugates(GEN T, long pinit)
    2311                 :            : {
    2312                 :        455 :   pari_sp av = avma;
    2313                 :        455 :   long c, nbtest, nbmax, n = degpol(T);
    2314                 :            :   ulong p;
    2315                 :            :   forprime_t S;
    2316                 :            : 
    2317         [ -  + ]:        455 :   if (n == 1) return 1;
    2318         [ +  + ]:        455 :   nbmax = (n < 10)? 20: (n<<1) + 1;
    2319                 :        455 :   nbtest = 0;
    2320                 :            : #if 0
    2321                 :            :   c = sturm(T); c = ugcd(c, n - c); /* too costly: finite primes are cheaper */
    2322                 :            : #else
    2323                 :        455 :   c = n;
    2324                 :            : #endif
    2325                 :        455 :   u_forprime_init(&S, pinit, ULONG_MAX);
    2326         [ +  - ]:       3759 :   while((p = u_forprime_next(&S)))
    2327                 :            :   {
    2328                 :       3759 :     GEN L, Tp = ZX_to_Flx(T,p);
    2329                 :            :     long i, nb;
    2330         [ +  + ]:       3759 :     if (!Flx_is_squarefree(Tp, p)) continue;
    2331                 :            :     /* unramified */
    2332                 :       3080 :     nbtest++;
    2333                 :       3080 :     L = Flx_nbfact_by_degree(Tp, &nb, p); /* L[i] = #factors of degree i */
    2334         [ +  + ]:       3080 :     if (L[n/nb] == nb) {
    2335 [ +  + ][ +  + ]:       1995 :       if (c == n && nbtest > 10) break; /* probably Galois */
    2336                 :            :     }
    2337                 :            :     else
    2338                 :            :     {
    2339                 :       1085 :       c = ugcd(c, L[1]);
    2340         [ +  + ]:       7476 :       for (i = 2; i <= n; i++)
    2341 [ +  + ][ +  + ]:       6727 :         if (L[i]) { c = ugcd(c, L[i]*i); if (c == 1) break; }
    2342         [ +  + ]:       1085 :       if (c == 1) break;
    2343                 :            :     }
    2344         [ +  + ]:       2737 :     if (nbtest == nbmax) break;
    2345         [ -  + ]:       2625 :     if (DEBUGLEVEL >= 6)
    2346                 :          0 :       err_printf("NumberOfConjugates [%ld]:c=%ld,p=%ld\n", nbtest,c,p);
    2347                 :       3080 :     avma = av;
    2348                 :            :   }
    2349         [ -  + ]:        455 :   if (DEBUGLEVEL >= 2) err_printf("NumberOfConjugates:c=%ld,p=%ld\n", c, p);
    2350                 :        455 :   avma = av; return c;
    2351                 :            : }
    2352                 :            : static GEN
    2353                 :          0 : galoisconj4(GEN nf, GEN d)
    2354                 :            : {
    2355                 :          0 :   pari_sp av = avma;
    2356                 :            :   GEN G, T;
    2357                 :          0 :   G = galoisconj4_main(nf, d, 0);
    2358         [ #  # ]:          0 :   if (typ(G) != t_INT) return G; /* Success */
    2359                 :          0 :   avma = av; T = get_nfpol(nf, &nf);
    2360                 :          0 :   G = cgetg(2, t_COL); gel(G,1) = pol_x(varn(T)); return G; /* Fail */
    2361                 :            : 
    2362                 :            : }
    2363                 :            : 
    2364                 :            : /* d multipllicative bound for the automorphism's denominators */
    2365                 :            : GEN
    2366                 :       1372 : galoisconj(GEN nf, GEN d)
    2367                 :            : {
    2368                 :       1372 :   pari_sp av = avma;
    2369                 :       1372 :   GEN G = galoisconj4_main(nf, d, 0);
    2370         [ +  + ]:       1372 :   if (typ(G) != t_INT) return G; /* Success */
    2371                 :       1372 :   avma = av; return galoisconj1(nf);
    2372                 :            : }
    2373                 :            : 
    2374                 :            : GEN
    2375                 :         42 : galoisconj0(GEN nf, long flag, GEN d, long prec)
    2376                 :            : {
    2377   [ +  +  -  -  :         42 :   switch(flag) {
                      - ]
    2378                 :         35 :     case 0: return galoisconj(nf, d);
    2379                 :          7 :     case 1: return galoisconj1(nf);
    2380                 :          0 :     case 2: return galoisconj2(nf, prec);
    2381                 :          0 :     case 4: return galoisconj4(nf, d);
    2382                 :            :   }
    2383                 :          0 :   pari_err_FLAG("nfgaloisconj");
    2384                 :         42 :   return NULL; /*not reached*/
    2385                 :            : }
    2386                 :            : 
    2387                 :            : /******************************************************************************/
    2388                 :            : /* Galois theory related algorithms                                           */
    2389                 :            : /******************************************************************************/
    2390                 :            : GEN
    2391                 :       2023 : checkgal(GEN gal)
    2392                 :            : {
    2393         [ -  + ]:       2023 :   if (typ(gal) == t_POL) pari_err_TYPE("checkgal [apply galoisinit first]",gal);
    2394 [ +  + ][ -  + ]:       2023 :   if (typ(gal) != t_VEC || lg(gal) != 9) pari_err_TYPE("checkgal",gal);
    2395                 :       2016 :   return gal;
    2396                 :            : }
    2397                 :            : 
    2398                 :            : GEN
    2399                 :        602 : galoisinit(GEN nf, GEN den)
    2400                 :            : {
    2401                 :        602 :   GEN G = galoisconj4_main(nf, den, 1);
    2402         [ +  + ]:        588 :   return (typ(G) == t_INT)? gen_0: G;
    2403                 :            : }
    2404                 :            : 
    2405                 :            : static GEN
    2406                 :       1225 : galoispermtopol_i(GEN gal, GEN perm, GEN mod, GEN mod2)
    2407                 :            : {
    2408                 :       1225 :   long t = typ(perm), i;
    2409                 :            :   GEN v;
    2410      [ +  +  - ]:       1225 :   switch (t)
    2411                 :            :   {
    2412                 :            :   case t_VECSMALL:
    2413                 :            :   {
    2414                 :        994 :     long v = varn(gel(gal,1));
    2415                 :        994 :     return permtopol(perm, gal_get_roots(gal), gal_get_invvdm(gal),
    2416                 :            :                            gal_get_den(gal), mod, mod2, v);
    2417                 :            :   }
    2418                 :            :   case t_VEC: case t_COL: case t_MAT:
    2419                 :        231 :     v = cgetg(lg(perm), t);
    2420         [ -  + ]:        231 :     if (DEBUGLEVEL>=4) err_printf("GaloisPermToPol:");
    2421         [ +  + ]:        693 :     for (i = 1; i < lg(v); i++)
    2422                 :            :     {
    2423                 :        462 :       gel(v,i) = galoispermtopol_i(gal, gel(perm,i), mod, mod2);
    2424         [ -  + ]:        462 :       if (DEBUGLEVEL>=4) err_printf("%ld ",i);
    2425                 :            :     }
    2426         [ -  + ]:        231 :     if (DEBUGLEVEL>=4) err_printf("\n");
    2427                 :        231 :     return v;
    2428                 :            :   }
    2429                 :          0 :   pari_err_TYPE("galoispermtopol", perm);
    2430                 :       1225 :   return NULL; /* not reached */
    2431                 :            : }
    2432                 :            : 
    2433                 :            : GEN
    2434                 :        763 : galoispermtopol(GEN gal, GEN perm)
    2435                 :            : {
    2436                 :        763 :   pari_sp av = avma;
    2437                 :            :   GEN mod, mod2;
    2438                 :        763 :   gal = checkgal(gal);
    2439                 :        763 :   mod = gal_get_mod(gal);
    2440                 :        763 :   mod2 = shifti(mod,-1);
    2441                 :        763 :   return gerepilecopy(av, galoispermtopol_i(gal, perm, mod, mod2));
    2442                 :            : }
    2443                 :            : 
    2444                 :            : GEN
    2445                 :         42 : galoiscosets(GEN O, GEN perm)
    2446                 :            : {
    2447                 :         42 :   long i, j, k, u, f, l = lg(O);
    2448                 :         42 :   GEN RC, C = cgetg(l,t_VECSMALL), o = gel(O,1);
    2449                 :         42 :   pari_sp av = avma;
    2450                 :         42 :   f = lg(o); u = o[1]; RC = zero_zv(lg(perm)-1);
    2451         [ +  + ]:        175 :   for(i=1,j=1; j<l; i++)
    2452                 :            :   {
    2453                 :        133 :     GEN p = gel(perm,i);
    2454         [ +  + ]:        133 :     if (RC[ p[u] ]) continue;
    2455         [ +  + ]:        315 :     for(k=1; k<f; k++) RC[ p[ o[k] ] ] = 1;
    2456                 :         91 :     C[j++] = i;
    2457                 :            :   }
    2458                 :         42 :   avma = av; return C;
    2459                 :            : }
    2460                 :            : 
    2461                 :            : static GEN
    2462                 :         42 : fixedfieldfactor(GEN L, GEN O, GEN perm, GEN M, GEN den, GEN mod, GEN mod2,
    2463                 :            :                  long x,long y)
    2464                 :            : {
    2465                 :         42 :   pari_sp ltop = avma;
    2466                 :         42 :   long i, j, k, l = lg(O), lo = lg(gel(O,1));
    2467                 :         42 :   GEN V, res, cosets = galoiscosets(O,perm), F = cgetg(lo+1,t_COL);
    2468                 :            : 
    2469                 :         42 :   gel(F, lo) = gen_1;
    2470         [ -  + ]:         42 :   if (DEBUGLEVEL>=4) err_printf("GaloisFixedField:cosets=%Ps \n",cosets);
    2471         [ -  + ]:         42 :   if (DEBUGLEVEL>=6) err_printf("GaloisFixedField:den=%Ps mod=%Ps \n",den,mod);
    2472                 :         42 :   V = cgetg(l,t_COL); res = cgetg(l,t_VEC);
    2473         [ +  + ]:        133 :   for (i = 1; i < l; i++)
    2474                 :            :   {
    2475                 :         91 :     pari_sp av = avma;
    2476                 :         91 :     GEN G = cgetg(l,t_VEC), Lp = vecpermute(L, gel(perm, cosets[i]));
    2477         [ +  + ]:        322 :     for (k = 1; k < l; k++)
    2478                 :        231 :       gel(G,k) = FpV_roots_to_pol(vecpermute(Lp, gel(O,k)), mod, x);
    2479         [ +  + ]:        315 :     for (j = 1; j < lo; j++)
    2480                 :            :     {
    2481         [ +  + ]:        700 :       for(k = 1; k < l; k++) gel(V,k) = gmael(G,k,j+1);
    2482                 :        224 :       gel(F,j) = vectopol(V, M, den, mod, mod2, y);
    2483                 :            :     }
    2484                 :         91 :     gel(res,i) = gerepileupto(av,gtopolyrev(F,x));
    2485                 :            :   }
    2486                 :         42 :   return gerepileupto(ltop,res);
    2487                 :            : }
    2488                 :            : 
    2489                 :            : static void
    2490                 :        994 : chk_perm(GEN perm, long n)
    2491                 :            : {
    2492 [ +  - ][ -  + ]:        994 :   if (typ(perm) != t_VECSMALL || lg(perm)!=n+1)
    2493                 :          0 :     pari_err_TYPE("galoisfixedfield", perm);
    2494                 :        994 : }
    2495                 :            : 
    2496                 :            : static int
    2497                 :       1281 : is_group(GEN g)
    2498                 :            : {
    2499 [ +  + ][ +  + ]:       2555 :   return typ(g)==t_VEC && lg(g)==3 && typ(gel(g,1))==t_VEC
    2500 [ +  + ][ +  - ]:       2555 :       && typ(gel(g,2))==t_VECSMALL;
    2501                 :            : }
    2502                 :            : 
    2503                 :            : GEN
    2504                 :        749 : galoisfixedfield(GEN gal, GEN perm, long flag, long y)
    2505                 :            : {
    2506                 :        749 :   pari_sp ltop = avma;
    2507                 :            :   GEN T, L, P, S, PL, O, res, mod, mod2;
    2508                 :            :   long vT, n, i;
    2509 [ +  - ][ -  + ]:        749 :   if (flag<0 || flag>2) pari_err_FLAG("galoisfixedfield");
    2510                 :        749 :   gal = checkgal(gal); T = gal_get_pol(gal);
    2511                 :        749 :   vT = varn(T);
    2512                 :        749 :   L = gal_get_roots(gal); n = lg(L)-1;
    2513                 :        749 :   mod = gal_get_mod(gal);
    2514         [ +  + ]:        749 :   if (typ(perm) == t_VEC)
    2515                 :            :   {
    2516         [ -  + ]:        525 :     if (is_group(perm)) perm = gel(perm, 1);
    2517         [ +  + ]:       1295 :     for (i = 1; i < lg(perm); i++) chk_perm(gel(perm,i), n);
    2518                 :        525 :     O = vecperm_orbits(perm, n);
    2519                 :            :   }
    2520                 :            :   else
    2521                 :            :   {
    2522                 :        224 :     chk_perm(perm, n);
    2523                 :        224 :     O = perm_cycles(perm);
    2524                 :            :   }
    2525                 :            : 
    2526                 :            :   {
    2527                 :        749 :     GEN OL= fixedfieldorbits(O,L);
    2528                 :        749 :     GEN V = fixedfieldsympol(OL, mod, gal_get_p(gal), NULL, vT);
    2529                 :        749 :     PL= gel(V,2);
    2530                 :        749 :     P = gel(V,3);
    2531                 :            :   }
    2532         [ +  + ]:        749 :   if (flag==1) return gerepileupto(ltop,P);
    2533                 :        539 :   mod2 = shifti(mod,-1);
    2534                 :        539 :   S = fixedfieldinclusion(O, PL);
    2535                 :        539 :   S = vectopol(S, gal_get_invvdm(gal), gal_get_den(gal), mod, mod2, vT);
    2536         [ +  + ]:        539 :   if (flag==0)
    2537                 :        497 :     res = cgetg(3, t_VEC);
    2538                 :            :   else
    2539                 :            :   {
    2540                 :            :     GEN PM, Pden;
    2541                 :            :     struct galois_borne Pgb;
    2542                 :         42 :     long val = itos(gal_get_e(gal));
    2543                 :         42 :     Pgb.l = gal_get_p(gal);
    2544                 :         42 :     Pden = galoisborne(P, NULL, &Pgb, degpol(T)/degpol(P));
    2545         [ -  + ]:         42 :     if (Pgb.valabs > val)
    2546                 :            :     {
    2547         [ #  # ]:          0 :       if (DEBUGLEVEL>=4)
    2548                 :          0 :         err_printf("GaloisConj:increase p-adic prec by %ld.\n", Pgb.valabs-val);
    2549                 :          0 :       PL = ZpX_liftroots(P, PL, Pgb.l, Pgb.valabs);
    2550                 :          0 :       L  = ZpX_liftroots(T, L, Pgb.l, Pgb.valabs);
    2551                 :          0 :       mod = Pgb.ladicabs; mod2 = shifti(mod,-1);
    2552                 :            :     }
    2553                 :         42 :     PM = vandermondeinversemod(PL, P, Pden, mod);
    2554         [ +  + ]:         42 :     if (y < 0) y = fetch_user_var("y");
    2555         [ -  + ]:         42 :     if (varncmp(y, vT) <= 0)
    2556                 :          0 :       pari_err_PRIORITY("galoisfixedfield", T, "<=", y);
    2557                 :         42 :     res = cgetg(4, t_VEC);
    2558                 :         42 :     gel(res,3) = fixedfieldfactor(L,O,gal_get_group(gal), PM,Pden,mod,mod2,vT,y);
    2559                 :            :   }
    2560                 :        539 :   gel(res,1) = gcopy(P);
    2561                 :        539 :   gel(res,2) = gmodulo(S, T);
    2562                 :        749 :   return gerepileupto(ltop, res);
    2563                 :            : }
    2564                 :            : 
    2565                 :            : /* gal a galois group output the underlying wss group */
    2566                 :            : GEN
    2567                 :        406 : galois_group(GEN gal) { return mkvec2(gal_get_gen(gal), gal_get_orders(gal)); }
    2568                 :            : 
    2569                 :            : GEN
    2570                 :        756 : checkgroup(GEN g, GEN *S)
    2571                 :            : {
    2572         [ +  + ]:        756 :   if (is_group(g)) { *S = NULL; return g; }
    2573                 :        399 :   g  = checkgal(g);
    2574                 :        749 :   *S = gal_get_group(g); return galois_group(g);
    2575                 :            : }
    2576                 :            : 
    2577                 :            : GEN
    2578                 :        168 : galoisisabelian(GEN gal, long flag)
    2579                 :            : {
    2580                 :        168 :   pari_sp av = avma;
    2581                 :        168 :   GEN S, G = checkgroup(gal,&S);
    2582         [ +  + ]:        168 :   if (!group_isabelian(G)) { avma=av; return gen_0; }
    2583   [ +  +  +  - ]:        147 :   switch(flag)
    2584                 :            :   {
    2585                 :         49 :     case 0: return gerepileupto(av, group_abelianHNF(G,S));
    2586                 :         49 :     case 1: avma=av; return gen_1;
    2587                 :         49 :     case 2: return gerepileupto(av, group_abelianSNF(G,S));
    2588                 :          0 :     default: pari_err_FLAG("galoisisabelian");
    2589                 :            :   }
    2590                 :        168 :   return NULL; /* not reached */
    2591                 :            : }
    2592                 :            : 
    2593                 :            : long
    2594                 :         56 : galoisisnormal(GEN gal, GEN sub)
    2595                 :            : {
    2596                 :         56 :   pari_sp av = avma;
    2597                 :         56 :   GEN S, G = checkgroup(gal, &S), H = checkgroup(sub, &S);
    2598                 :         56 :   long res = group_subgroup_isnormal(G, H);
    2599                 :         56 :   avma = av;
    2600                 :         56 :   return res;
    2601                 :            : }
    2602                 :            : 
    2603                 :            : GEN
    2604                 :         56 : galoissubgroups(GEN gal)
    2605                 :            : {
    2606                 :         56 :   pari_sp av = avma;
    2607                 :         56 :   GEN S, G = checkgroup(gal,&S);
    2608                 :         56 :   return gerepileupto(av, group_subgroups(G));
    2609                 :            : }
    2610                 :            : 
    2611                 :            : GEN
    2612                 :         42 : galoissubfields(GEN G, long flag, long v)
    2613                 :            : {
    2614                 :         42 :   pari_sp av = avma;
    2615                 :         42 :   GEN L = galoissubgroups(G);
    2616                 :         42 :   long i, l = lg(L);
    2617                 :         42 :   GEN S = cgetg(l, t_VEC);
    2618         [ +  + ]:        567 :   for (i = 1; i < l; ++i) gel(S,i) = galoisfixedfield(G, gmael(L,i,1), flag, v);
    2619                 :         42 :   return gerepileupto(av, S);
    2620                 :            : }
    2621                 :            : 
    2622                 :            : GEN
    2623                 :         28 : galoisexport(GEN gal, long format)
    2624                 :            : {
    2625                 :         28 :   pari_sp av = avma;
    2626                 :         28 :   GEN S, G = checkgroup(gal,&S);
    2627                 :         28 :   return gerepileupto(av, group_export(G,format));
    2628                 :            : }
    2629                 :            : 
    2630                 :            : GEN
    2631                 :        392 : galoisidentify(GEN gal)
    2632                 :            : {
    2633                 :        392 :   pari_sp av = avma;
    2634                 :        392 :   GEN S, G = checkgroup(gal,&S);
    2635                 :        385 :   long idx = group_ident(G,S), card = group_order(G);
    2636                 :        385 :   avma = av; return mkvec2s(card, idx);
    2637                 :            : }

Generated by: LCOV version 1.9