Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - galconj.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25406-bf255ab81b) Lines: 1722 2132 80.8 %
Date: 2020-06-04 05:59:24 Functions: 111 145 76.6 %
Legend: Lines: hit not hit

          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             : static int
      23         280 : is2sparse(GEN x)
      24             : {
      25         280 :   long i, l = lg(x);
      26         280 :   if (odd(l-3)) return 0;
      27         826 :   for(i=3; i<l; i+=2)
      28         651 :     if (signe(gel(x,i))) return 0;
      29         175 :   return 1;
      30             : }
      31             : 
      32             : static GEN
      33         917 : galoisconj1(GEN nf)
      34             : {
      35         917 :   GEN x = get_nfpol(nf, &nf), f = nf? nf : x, y, z;
      36         917 :   long i, lz, v = varn(x), nbmax;
      37         917 :   pari_sp av = avma;
      38         917 :   RgX_check_ZX(x, "nfgaloisconj");
      39         917 :   nbmax = numberofconjugates(x, 2);
      40         917 :   if (nbmax==1) retmkcol(pol_x(v));
      41         378 :   if (nbmax==2 && is2sparse(x))
      42             :   {
      43         175 :     GEN res = cgetg(3,t_COL);
      44         175 :     gel(res,1) = deg1pol_shallow(gen_m1, gen_0, v);
      45         175 :     gel(res,2) = pol_x(v);
      46         175 :     return res;
      47             :   }
      48         203 :   x = leafcopy(x); setvarn(x, fetch_var_higher());
      49         203 :   z = nfroots(f, x); lz = lg(z);
      50         203 :   y = cgetg(lz, t_COL);
      51        1008 :   for (i = 1; i < lz; i++)
      52             :   {
      53         805 :     GEN t = lift(gel(z,i));
      54         805 :     if (typ(t) == t_POL) setvarn(t, v);
      55         805 :     gel(y,i) = t;
      56             :   }
      57         203 :   (void)delete_var();
      58         203 :   return gerepileupto(av, y);
      59             : }
      60             : 
      61             : /*************************************************************************/
      62             : /**                                                                     **/
      63             : /**                           GALOISCONJ4                               **/
      64             : /**                                                                     **/
      65             : /*************************************************************************/
      66             : /*DEBUGLEVEL:
      67             :   1: timing
      68             :   2: outline
      69             :   4: complete outline
      70             :   6: detail
      71             :   7: memory
      72             :   9: complete detail
      73             : */
      74             : struct galois_borne {
      75             :   GEN  l;
      76             :   long valsol;
      77             :   long valabs;
      78             :   GEN  bornesol;
      79             :   GEN  ladicsol;
      80             :   GEN  ladicabs;
      81             :   GEN  dis;
      82             : };
      83             : struct galois_lift {
      84             :   GEN  T;
      85             :   GEN  den;
      86             :   GEN  p;
      87             :   GEN  L;
      88             :   GEN  Lden;
      89             :   long e;
      90             :   GEN  Q;
      91             :   GEN  TQ;
      92             :   struct galois_borne *gb;
      93             : };
      94             : struct galois_testlift {
      95             :   long n;
      96             :   long f;
      97             :   long g;
      98             :   GEN  bezoutcoeff;
      99             :   GEN  pauto;
     100             :   GEN  C;
     101             :   GEN  Cd;
     102             : };
     103             : struct galois_test { /* data for permutation test */
     104             :   GEN order; /* order of tests pour galois_test_perm */
     105             :   GEN borne, lborne; /* coefficient bounds */
     106             :   GEN ladic;
     107             :   GEN PV; /* NULL or vector of test matrices (Vmatrix) */
     108             :   GEN TM; /* transpose of M */
     109             :   GEN L; /* p-adic roots, known mod ladic */
     110             :   GEN M; /* vandermonde inverse */
     111             : };
     112             : /* result of the study of Frobenius degrees */
     113             : enum ga_code {ga_all_normal=1,ga_ext_2=2,ga_non_wss=4,
     114             :               ga_all_nilpotent=8,ga_easy=16};
     115             : struct galois_analysis {
     116             :   long p; /* prime to be lifted */
     117             :   long deg; /* degree of the lift */
     118             :   long ord;
     119             :   long l; /* l: prime number such that T is totally split mod l */
     120             :   long p4;
     121             :   long group;
     122             : };
     123             : struct galois_frobenius {
     124             :   long p;
     125             :   long fp;
     126             :   long deg;
     127             :   GEN Tmod;
     128             :   GEN psi;
     129             : };
     130             : 
     131             : /* given complex roots L[i], i <= n of some monic T in C[X], return
     132             :  * the T'(L[i]), computed stably via products of differences */
     133             : static GEN
     134       10587 : vandermondeinverseprep(GEN L)
     135             : {
     136       10587 :   long i, j, n = lg(L);
     137             :   GEN V;
     138       10587 :   V = cgetg(n, t_VEC);
     139       66814 :   for (i = 1; i < n; i++)
     140             :   {
     141       56227 :     pari_sp ltop = avma;
     142       56227 :     GEN W = cgetg(n-1,t_VEC);
     143       56227 :     long k = 1;
     144      991206 :     for (j = 1; j < n; j++)
     145      934979 :       if (i != j) gel(W, k++) = gsub(gel(L,i),gel(L,j));
     146       56227 :     gel(V,i) = gerepileupto(ltop, RgV_prod(W));
     147             :   }
     148       10587 :   return V;
     149             : }
     150             : 
     151             : /* Compute the inverse of the van der Monde matrix of T multiplied by den */
     152             : GEN
     153       10482 : vandermondeinverse(GEN L, GEN T, GEN den, GEN prep)
     154             : {
     155       10482 :   pari_sp ltop = avma;
     156       10482 :   long i, n = lg(L)-1;
     157             :   GEN M, P;
     158       10482 :   if (!prep) prep = vandermondeinverseprep(L);
     159       10482 :   if (den && !equali1(den)) T = RgX_Rg_mul(T,den);
     160       10482 :   M = cgetg(n+1, t_MAT);
     161       64917 :   for (i = 1; i <= n; i++)
     162             :   {
     163       54435 :     P = RgX_Rg_div(RgX_div_by_X_x(T, gel(L,i), NULL), gel(prep,i));
     164       54435 :     gel(M,i) = RgX_to_RgC(P,n);
     165             :   }
     166       10482 :   return gerepilecopy(ltop, M);
     167             : }
     168             : 
     169             : /* #r = r1 + r2 */
     170             : GEN
     171          78 : embed_roots(GEN ro, long r1)
     172             : {
     173          78 :   long r2 = lg(ro)-1-r1;
     174             :   GEN L;
     175          78 :   if (!r2) L = ro;
     176             :   else
     177             :   {
     178          50 :     long i,j, N = r1+2*r2;
     179          50 :     L = cgetg(N+1, t_VEC);
     180         159 :     for (i = 1; i <= r1; i++) gel(L,i) = gel(ro,i);
     181         410 :     for (j = i; j <= N; i++)
     182             :     {
     183         360 :       GEN z = gel(ro,i);
     184         360 :       gel(L,j++) = z;
     185         360 :       gel(L,j++) = mkcomplex(gel(z,1), gneg(gel(z,2)));
     186             :     }
     187             :   }
     188          78 :   return L;
     189             : }
     190             : GEN
     191       35399 : embed_disc(GEN z, long r1, long prec)
     192             : {
     193       35399 :   pari_sp av = avma;
     194       35399 :   GEN t = real_1(prec);
     195       35399 :   long i, j, n = lg(z)-1, r2 = n-r1;
     196      130431 :   for (i = 1; i < r1; i++)
     197             :   {
     198       95032 :     GEN zi = gel(z,i);
     199      665042 :     for (j = i+1; j <= r1; j++) t = gmul(t, gsub(zi, gel(z,j)));
     200             :   }
     201      193242 :   for (j = r1+1; j <= n; j++)
     202             :   {
     203      157843 :     GEN zj = gel(z,j), a = gel(zj,1), b = gel(zj,2), b2 = gsqr(b);
     204      174601 :     for (i = 1; i <= r1; i++)
     205             :     {
     206       16758 :       GEN zi = gel(z,i);
     207       16758 :       t = gmul(t, gadd(gsqr(gsub(zi, a)), b2));
     208             :     }
     209      157843 :     t = gmul(t, b);
     210             :   }
     211       35399 :   if (r2) t = gmul2n(t, r2);
     212       35399 :   if (r2 > 1)
     213             :   {
     214       24164 :     GEN T = real_1(prec);
     215      156975 :     for (i = r1+1; i < n; i++)
     216             :     {
     217      132811 :       GEN zi = gel(z,i), a = gel(zi,1), b = gel(zi,2);
     218      864591 :       for (j = i+1; j <= n; j++)
     219             :       {
     220      731780 :         GEN zj = gel(z,j), c = gel(zj,1), d = gel(zj,2);
     221      731780 :         GEN f = gsqr(gsub(a,c)), g = gsqr(gsub(b,d)), h = gsqr(gadd(b,d));
     222      731780 :         T = gmul(T, gmul(gadd(f,g), gadd(f,h)));
     223             :       }
     224             :     }
     225       24164 :     t = gmul(t, T);
     226             :   }
     227       35399 :   t = gsqr(t);
     228       35399 :   if (odd(r2)) t = gneg(t);
     229       35399 :   return gerepileupto(av, t);
     230             : }
     231             : 
     232             : /* Compute bound for the coefficients of automorphisms.
     233             :  * T a ZX, den a t_INT denominator or NULL */
     234             : GEN
     235       10587 : initgaloisborne(GEN T, GEN den, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis)
     236             : {
     237             :   GEN L, prep, nf, r;
     238             :   pari_timer ti;
     239             : 
     240       10587 :   if (DEBUGLEVEL>=4) timer_start(&ti);
     241       10587 :   T = get_nfpol(T, &nf);
     242       10587 :   r = nf ? nf_get_roots(nf) : NULL;
     243       10587 :   if (nf &&  precision(gel(r, 1)) >= prec)
     244          78 :     L = embed_roots(r, nf_get_r1(nf));
     245             :   else
     246       10509 :     L = QX_complex_roots(T, prec);
     247       10587 :   if (DEBUGLEVEL>=4) timer_printf(&ti,"roots");
     248       10587 :   prep = vandermondeinverseprep(L);
     249       10587 :   if (!den || ptdis)
     250             :   {
     251        3346 :     GEN res = RgV_prod(gabs(prep,prec));
     252        3346 :     GEN D = ZX_disc_all(T, 1 + expi(ceil_safe(res))); /* +1 if inaccurate res */
     253        3346 :     if (ptdis) *ptdis = D;
     254        3346 :     if (!den) den = indexpartial(T,D);
     255             :   }
     256       10587 :   if (ptprep) *ptprep = prep;
     257       10587 :   *ptL = L; return den;
     258             : }
     259             : 
     260             : /* ||| M ||| with respect to || x ||_oo, M t_MAT */
     261             : GEN
     262       26714 : matrixnorm(GEN M, long prec)
     263             : {
     264       26714 :   long i,j,m, l = lg(M);
     265       26714 :   GEN B = real_0(prec);
     266             : 
     267       26714 :   if (l == 1) return B;
     268       26714 :   m = lgcols(M);
     269       84989 :   for (i = 1; i < m; i++)
     270             :   {
     271       58275 :     GEN z = gabs(gcoeff(M,i,1), prec);
     272      772104 :     for (j = 2; j < l; j++) z = gadd(z, gabs(gcoeff(M,i,j), prec));
     273       58275 :     if (gcmp(z, B) > 0) B = z;
     274             :   }
     275       26714 :   return B;
     276             : }
     277             : 
     278             : static GEN
     279        3241 : galoisborne(GEN T, GEN dn, struct galois_borne *gb, long d)
     280             : {
     281             :   pari_sp ltop, av2;
     282             :   GEN borne, borneroots, bornetrace, borneabs;
     283             :   long prec;
     284             :   GEN L, M, prep, den;
     285             :   pari_timer ti;
     286        3241 :   const long step=3;
     287             : 
     288        3241 :   prec = nbits2prec(bit_accuracy(ZX_max_lg(T)));
     289        3241 :   den = initgaloisborne(T,dn,prec, &L,&prep,&gb->dis);
     290        3241 :   if (!dn) dn = den;
     291        3241 :   ltop = avma;
     292        3241 :   if (DEBUGLEVEL>=4) timer_start(&ti);
     293        3241 :   M = vandermondeinverse(L, RgX_gtofp(T, prec), den, prep);
     294        3241 :   if (DEBUGLEVEL>=4) timer_printf(&ti,"vandermondeinverse");
     295        3241 :   borne = matrixnorm(M, prec);
     296        3241 :   borneroots = gsupnorm(L, prec); /*t_REAL*/
     297        3241 :   bornetrace = gmulsg((2*step)*degpol(T)/d,
     298        3241 :                       powru(borneroots, minss(degpol(T), step)));
     299        3241 :   borneroots = ceil_safe(gmul(borne, borneroots));
     300        3241 :   borneabs = ceil_safe(gmax_shallow(gmul(borne, bornetrace),
     301             :                                     powru(bornetrace, d)));
     302        3241 :   av2 = avma;
     303             :   /*We use d-1 test, so we must overlift to 2^BITS_IN_LONG*/
     304        3241 :   gb->valsol = logint(shifti(borneroots,2+BITS_IN_LONG), gb->l) + 1;
     305        3241 :   gb->valabs = logint(shifti(borneabs,2), gb->l) + 1;
     306        3241 :   gb->valabs = maxss(gb->valsol, gb->valabs);
     307        3241 :   if (DEBUGLEVEL >= 4)
     308           0 :     err_printf("GaloisConj: val1=%ld val2=%ld\n", gb->valsol, gb->valabs);
     309        3241 :   set_avma(av2);
     310        3241 :   gb->bornesol = gerepileuptoint(ltop, shifti(borneroots,1));
     311        3241 :   if (DEBUGLEVEL >= 9)
     312           0 :     err_printf("GaloisConj: Bound %Ps\n",borneroots);
     313        3241 :   gb->ladicsol = powiu(gb->l, gb->valsol);
     314        3241 :   gb->ladicabs = powiu(gb->l, gb->valabs);
     315        3241 :   return dn;
     316             : }
     317             : 
     318             : static GEN
     319        3059 : makeLden(GEN L,GEN den, struct galois_borne *gb)
     320        3059 : { return FpC_Fp_mul(L, den, gb->ladicsol); }
     321             : 
     322             : /* Initialize the galois_lift structure */
     323             : static void
     324        3122 : initlift(GEN T, GEN den, ulong p, GEN L, GEN Lden, struct galois_borne *gb, struct galois_lift *gl)
     325             : {
     326             :   pari_sp av;
     327             :   long e;
     328        3122 :   gl->gb = gb;
     329        3122 :   gl->T = T;
     330        3122 :   gl->den = is_pm1(den)? gen_1: den;
     331        3122 :   gl->p = utoipos(p);
     332        3122 :   gl->L = L;
     333        3122 :   gl->Lden = Lden;
     334        3122 :   av = avma;
     335        3122 :   e = logint(shifti(gb->bornesol, 2+BITS_IN_LONG), gl->p) + 1;
     336        3122 :   set_avma(av);
     337        3122 :   if (e < 2) e = 2;
     338        3122 :   gl->e = e;
     339        3122 :   gl->Q = powuu(p, e);
     340        3122 :   gl->TQ = FpX_red(T,gl->Q);
     341        3122 : }
     342             : 
     343             : /* Check whether f is (with high probability) a solution and compute its
     344             :  * permutation */
     345             : static int
     346        7426 : poltopermtest(GEN f, struct galois_lift *gl, GEN pf)
     347             : {
     348             :   pari_sp av;
     349        7426 :   GEN fx, fp, B = gl->gb->bornesol;
     350             :   long i, j, ll;
     351       53352 :   for (i = 2; i < lg(f); i++)
     352       47859 :     if (abscmpii(gel(f,i),B) > 0)
     353             :     {
     354        1933 :       if (DEBUGLEVEL>=4) err_printf("GaloisConj: Solution too large.\n");
     355        1933 :       if (DEBUGLEVEL>=8) err_printf("f=%Ps\n borne=%Ps\n",f,B);
     356        1933 :       return 0;
     357             :     }
     358        5493 :   ll = lg(gl->L);
     359        5493 :   fp = const_vecsmall(ll-1, 1); /* left on stack */
     360        5493 :   av = avma;
     361       60718 :   for (i = 1; i < ll; i++, set_avma(av))
     362             :   {
     363       55245 :     fx = FpX_eval(f, gel(gl->L,i), gl->gb->ladicsol);
     364      743849 :     for (j = 1; j < ll; j++)
     365      743829 :       if (fp[j] && equalii(fx, gel(gl->Lden,j))) { pf[i]=j; fp[j]=0; break; }
     366       55245 :     if (j == ll) return 0;
     367             :   }
     368        5473 :   return 1;
     369             : }
     370             : 
     371             : static long
     372        6550 : galoisfrobeniustest(GEN aut, struct galois_lift *gl, GEN frob)
     373             : {
     374        6550 :   pari_sp av = avma;
     375        6550 :   GEN tlift = aut;
     376        6550 :   if (gl->den != gen_1) tlift = FpX_Fp_mul(tlift, gl->den, gl->Q);
     377        6550 :   tlift = FpX_center_i(tlift, gl->Q, shifti(gl->Q,-1));
     378        6550 :   return gc_long(av, poltopermtest(tlift, gl, frob));
     379             : }
     380             : 
     381             : static GEN
     382       13265 : monoratlift(void *E, GEN S, GEN q)
     383             : {
     384       13265 :   pari_sp ltop = avma;
     385       13265 :   struct galois_lift *gl = (struct galois_lift *) E;
     386       13265 :   GEN qm1 = sqrti(shifti(q,-2)), N = gl->Q;
     387       13265 :   GEN tlift = FpX_ratlift(S, q, qm1, qm1, gl->den);
     388       13265 :   if (tlift)
     389             :   {
     390        2154 :     pari_sp ltop = avma;
     391        2154 :     GEN frob = cgetg(lg(gl->L), t_VECSMALL);
     392        2154 :     if(DEBUGLEVEL>=4)
     393           0 :       err_printf("MonomorphismLift: trying early solution %Ps\n",tlift);
     394        2154 :     if (gl->den != gen_1)
     395        1594 :       tlift = FpX_Fp_mul(FpX_red(Q_muli_to_int(tlift, gl->den), N),
     396             :                          Fp_inv(gl->den, N), N);
     397        2154 :     if (galoisfrobeniustest(tlift, gl, frob))
     398             :     {
     399        2134 :       if(DEBUGLEVEL>=4) err_printf("MonomorphismLift: true early solution.\n");
     400        2134 :       return gerepilecopy(ltop, tlift);
     401             :     }
     402          20 :     if(DEBUGLEVEL>=4) err_printf("MonomorphismLift: false early solution.\n");
     403             :   }
     404       11131 :   set_avma(ltop);
     405       11131 :   return NULL;
     406             : }
     407             : 
     408             : static GEN
     409        3976 : monomorphismratlift(GEN P, GEN S, struct galois_lift *gl)
     410             : {
     411             :   pari_timer ti;
     412        3976 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
     413        3976 :   S = ZpX_ZpXQ_liftroot_ea(P,S,gl->T,gl->p, gl->e, (void*)gl, monoratlift);
     414        3976 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "monomorphismlift()");
     415        3976 :   return S;
     416             : }
     417             : 
     418             : /* Let T be a polynomial in Z[X] , p a prime number, S in Fp[X]/(T) so
     419             :  * that T(S)=0 [p,T]. Lift S in S_0 so that T(S_0)=0 [T,p^e]
     420             :  * Unclean stack */
     421             : static GEN
     422        3976 : automorphismlift(GEN S, struct galois_lift *gl)
     423             : {
     424        3976 :   return monomorphismratlift(gl->T, S, gl);
     425             : }
     426             : 
     427             : static GEN
     428        3122 : galoisdolift(struct galois_lift *gl)
     429             : {
     430        3122 :   pari_sp av = avma;
     431        3122 :   GEN Tp = FpX_red(gl->T, gl->p);
     432        3122 :   GEN S = FpX_Frobenius(Tp, gl->p);
     433        3122 :   return gerepileupto(av, automorphismlift(S, gl));
     434             : }
     435             : 
     436             : static GEN
     437         637 : galoisdoliftn(struct galois_lift *gl, long e)
     438             : {
     439         637 :   pari_sp av = avma;
     440         637 :   GEN Tp = FpX_red(gl->T, gl->p);
     441         637 :   GEN S = FpXQ_autpow(FpX_Frobenius(Tp, gl->p), e, Tp, gl->p);
     442         637 :   return gerepileupto(av, automorphismlift(S, gl));
     443             : }
     444             : 
     445             : static ulong
     446          70 : findpsi(GEN D, ulong pstart, GEN P, GEN S, long o, GEN *Tmod, GEN *Tpsi)
     447             : {
     448             :   forprime_t iter;
     449             :   ulong p;
     450          70 :   long n = degpol(P), i, j, g = n/o;
     451          70 :   GEN psi = cgetg(g+1, t_VECSMALL);
     452          70 :   u_forprime_init(&iter, pstart, ULONG_MAX);
     453        2191 :   while ((p = u_forprime_next(&iter)))
     454             :   {
     455             :     GEN F, Sp;
     456        2191 :     long gp = 0;
     457        2191 :     if (smodis(D, p) == 0)
     458         189 :       continue;
     459        2002 :     F = gel(Flx_factor(ZX_to_Flx(P, p), p), 1);
     460        2002 :     if (lg(F)-1 != g) continue;
     461         693 :     Sp = RgX_to_Flx(S, p);
     462        1778 :     for (j = 1; j <= g; j++)
     463             :     {
     464        1659 :       GEN Fj = gel(F, j);
     465        1659 :       GEN Sj = Flx_rem(Sp, Fj, p);
     466        1659 :       GEN A = Flxq_autpowers(Flx_Frobenius(Fj, p), o,  Fj, p);
     467        5670 :       for (i = 1; i <= o; i++)
     468        5096 :         if (gequal(Sj, gel(A,i+1)))
     469             :         {
     470        1085 :           psi[j] = i; break;
     471             :         }
     472        1659 :       if (i > o) break;
     473        1085 :       if (gp==0 && i==1) gp=j;
     474             :     }
     475         693 :     if (gp && j > g)
     476             :     {
     477             :       /* Normalize result so that psi[l]=1 */
     478          70 :       if (gp!=1)
     479             :       {
     480           7 :         swap(gel(F,1),gel(F,gp));
     481           7 :         lswap(uel(psi,1),uel(psi,gp));
     482             :       }
     483          70 :       *Tpsi = Flv_Fl_div(psi,psi[g],o);
     484          70 :       *Tmod = FlxV_to_ZXV(F);
     485          70 :       return p;
     486             :     }
     487             :   }
     488           0 :   return 0;
     489             : }
     490             : 
     491             : static void
     492         511 : inittestlift(GEN plift, GEN Tmod, struct galois_lift *gl,
     493             :              struct galois_testlift *gt)
     494             : {
     495             :   pari_timer ti;
     496         511 :   gt->n = lg(gl->L) - 1;
     497         511 :   gt->g = lg(Tmod) - 1;
     498         511 :   gt->f = gt->n / gt->g;
     499         511 :   gt->bezoutcoeff = bezout_lift_fact(gl->T, Tmod, gl->p, gl->e);
     500         511 :   if (DEBUGLEVEL >= 2) timer_start(&ti);
     501         511 :   gt->pauto = FpXQ_autpowers(plift, gt->f-1, gl->TQ, gl->Q);
     502         511 :   if (DEBUGLEVEL >= 2) timer_printf(&ti, "Frobenius power");
     503         511 : }
     504             : 
     505             : /* Explanation of the intheadlong technique:
     506             :  * Let C be a bound, B = BITS_IN_LONG, M > C*2^B a modulus and 0 <= a_i < M for
     507             :  * i=1,...,n where n < 2^B. We want to test if there exist k,l, |k| < C < M/2^B,
     508             :  * such that sum a_i = k + l*M
     509             :  * We write a_i*2^B/M = b_i+c_i with b_i integer and 0<=c_i<1, so that
     510             :  *   sum b_i - l*2^B = k*2^B/M - sum c_i
     511             :  * Since -1 < k*2^B/M < 1 and 0<=c_i<1, it follows that
     512             :  *   -n-1 < sum b_i - l*2^B < 1  i.e.  -n <= sum b_i -l*2^B <= 0
     513             :  * So we compute z = - sum b_i [mod 2^B] and check if 0 <= z <= n. */
     514             : 
     515             : /* Assume 0 <= x < mod. */
     516             : static ulong
     517     1064014 : intheadlong(GEN x, GEN mod)
     518             : {
     519     1064014 :   pari_sp av = avma;
     520     1064014 :   long res = (long) itou(divii(shifti(x,BITS_IN_LONG),mod));
     521     1064014 :   return gc_long(av,res);
     522             : }
     523             : static GEN
     524       32676 : vecheadlong(GEN W, GEN mod)
     525             : {
     526       32676 :   long i, l = lg(W);
     527       32676 :   GEN V = cgetg(l, t_VECSMALL);
     528     1066534 :   for(i=1; i<l; i++) V[i] = intheadlong(gel(W,i), mod);
     529       32676 :   return V;
     530             : }
     531             : static GEN
     532        2408 : matheadlong(GEN W, GEN mod)
     533             : {
     534        2408 :   long i, l = lg(W);
     535        2408 :   GEN V = cgetg(l,t_MAT);
     536       35084 :   for(i=1; i<l; i++) gel(V,i) = vecheadlong(gel(W,i), mod);
     537        2408 :   return V;
     538             : }
     539             : static ulong
     540       30156 : polheadlong(GEN P, long n, GEN mod)
     541             : {
     542       30156 :   return (lg(P)>n+2)? intheadlong(gel(P,n+2),mod): 0;
     543             : }
     544             : 
     545             : #define headlongisint(Z,N) (-(ulong)(Z)<=(ulong)(N))
     546             : 
     547             : static long
     548         504 : frobeniusliftall(GEN sg, long el, GEN *psi, struct galois_lift *gl,
     549             :                  struct galois_testlift *gt, GEN frob)
     550             : {
     551         504 :   pari_sp av, ltop2, ltop = avma;
     552         504 :   long i,j,k, c = lg(sg)-1, n = lg(gl->L)-1, m = gt->g, d = m / c;
     553             :   GEN pf, u, v, C, Cd, SG, cache;
     554         504 :   long N1, N2, R1, Ni, ord = gt->f, c_idx = gt->g-1;
     555             :   ulong headcache;
     556         504 :   long hop = 0;
     557             :   GEN NN, NQ;
     558             :   pari_timer ti;
     559             : 
     560         504 :   *psi = pf = cgetg(m, t_VECSMALL);
     561         504 :   ltop2 = avma;
     562         504 :   NN = diviiexact(mpfact(m), mului(c, powiu(mpfact(d), c)));
     563         504 :   if (DEBUGLEVEL >= 4)
     564           0 :     err_printf("GaloisConj: I will try %Ps permutations\n", NN);
     565         504 :   N1=10000000;
     566         504 :   NQ=divis_rem(NN,N1,&R1);
     567         504 :   if (abscmpiu(NQ,1000000000)>0)
     568             :   {
     569           0 :     pari_warn(warner,"Combinatorics too hard : would need %Ps tests!\n"
     570             :         "I will skip it, but it may induce an infinite loop",NN);
     571           0 :     *psi = NULL; return gc_long(ltop,0);
     572             :   }
     573         504 :   N2=itos(NQ); if(!N2) N1=R1;
     574         504 :   if (DEBUGLEVEL>=4) timer_start(&ti);
     575         504 :   set_avma(ltop2);
     576         504 :   C = gt->C;
     577         504 :   Cd= gt->Cd;
     578         504 :   v = FpXQ_mul(gel(gt->pauto, 1+el%ord), gel(gt->bezoutcoeff, m),gl->TQ,gl->Q);
     579         504 :   if (gl->den != gen_1) v = FpX_Fp_mul(v, gl->den, gl->Q);
     580         504 :   SG = cgetg(lg(sg),t_VECSMALL);
     581        1540 :   for(i=1; i<lg(SG); i++) SG[i] = (el*sg[i])%ord + 1;
     582         504 :   cache = cgetg(m+1,t_VECSMALL); cache[m] = polheadlong(v,1,gl->Q);
     583         504 :   headcache = polheadlong(v,2,gl->Q);
     584        1274 :   for (i = 1; i < m; i++) pf[i] = 1 + i/d;
     585         504 :   av = avma;
     586         504 :   for (Ni = 0, i = 0; ;i++)
     587             :   {
     588      292110 :     for (j = c_idx ; j > 0; j--)
     589             :     {
     590      230559 :       long h = SG[pf[j]];
     591      230559 :       if (!mael(C,h,j))
     592             :       {
     593        1050 :         pari_sp av3 = avma;
     594        1050 :         GEN r = FpXQ_mul(gel(gt->pauto,h), gel(gt->bezoutcoeff,j),gl->TQ,gl->Q);
     595        1050 :         if (gl->den != gen_1) r = FpX_Fp_mul(r, gl->den, gl->Q);
     596        1050 :         gmael(C,h,j) = gclone(r);
     597        1050 :         mael(Cd,h,j) = polheadlong(r,1,gl->Q);
     598        1050 :         set_avma(av3);
     599             :       }
     600      230559 :       uel(cache,j) = uel(cache,j+1)+umael(Cd,h,j);
     601             :     }
     602       61551 :     if (headlongisint(uel(cache,1),n))
     603             :     {
     604        2121 :       ulong head = headcache;
     605       30219 :       for (j = 1; j < m; j++) head += polheadlong(gmael(C,SG[pf[j]],j),2,gl->Q);
     606        2121 :       if (headlongisint(head,n))
     607             :       {
     608         483 :         u = v;
     609        1197 :         for (j = 1; j < m; j++) u = ZX_add(u, gmael(C,SG[pf[j]],j));
     610         483 :         u = FpX_center_i(FpX_red(u, gl->Q), gl->Q, shifti(gl->Q,-1));
     611         483 :         if (poltopermtest(u, gl, frob))
     612             :         {
     613         476 :           if (DEBUGLEVEL >= 4)
     614             :           {
     615           0 :             timer_printf(&ti, "");
     616           0 :             err_printf("GaloisConj: %d hops on %Ps tests\n",hop,addis(mulss(Ni,N1),i));
     617             :           }
     618         476 :           return gc_long(ltop2,1);
     619             :         }
     620           7 :         if (DEBUGLEVEL >= 4) err_printf("M");
     621             :       }
     622        1638 :       else hop++;
     623             :     }
     624       61075 :     if (DEBUGLEVEL >= 4 && i % maxss(N1/20, 1) == 0)
     625           0 :       timer_printf(&ti, "GaloisConj:Testing %Ps", addis(mulss(Ni,N1),i));
     626       61075 :     set_avma(av);
     627       61075 :     if (i == N1-1)
     628             :     {
     629          28 :       if (Ni==N2-1) N1 = R1;
     630          28 :       if (Ni==N2) break;
     631           0 :       Ni++; i = 0;
     632           0 :       if (DEBUGLEVEL>=4) timer_start(&ti);
     633             :     }
     634      168742 :     for (j = 2; j < m && pf[j-1] >= pf[j]; j++)
     635             :       /*empty*/; /* to kill clang Warning */
     636       99512 :     for (k = 1; k < j-k && pf[k] != pf[j-k]; k++) { lswap(pf[k], pf[j-k]); }
     637      107751 :     for (k = j - 1; pf[k] >= pf[j]; k--)
     638             :       /*empty*/;
     639       61047 :     lswap(pf[j], pf[k]); c_idx = j;
     640             :   }
     641          28 :   if (DEBUGLEVEL>=4) err_printf("GaloisConj: not found, %d hops \n",hop);
     642          28 :   *psi = NULL; return gc_long(ltop,0);
     643             : }
     644             : 
     645             : /* Compute the test matrix for the i-th line of V. Clone. */
     646             : static GEN
     647        2408 : Vmatrix(long i, struct galois_test *td)
     648             : {
     649        2408 :   pari_sp av = avma;
     650        2408 :   GEN m = gclone( matheadlong(FpC_FpV_mul(td->L, row(td->M,i), td->ladic), td->ladic));
     651        2408 :   set_avma(av); return m;
     652             : }
     653             : 
     654             : /* Initialize galois_test */
     655             : static void
     656        2184 : inittest(GEN L, GEN M, GEN borne, GEN ladic, struct galois_test *td)
     657             : {
     658        2184 :   long i, n = lg(L)-1;
     659        2184 :   GEN p = cgetg(n+1, t_VECSMALL);
     660        2184 :   if (DEBUGLEVEL >= 8) err_printf("GaloisConj: Init Test\n");
     661        2184 :   td->order = p;
     662       22694 :   for (i = 1; i <= n-2; i++) p[i] = i+2;
     663        2184 :   p[n-1] = 1; p[n] = 2;
     664        2184 :   td->borne = borne;
     665        2184 :   td->lborne = subii(ladic, borne);
     666        2184 :   td->ladic = ladic;
     667        2184 :   td->L = L;
     668        2184 :   td->M = M;
     669        2184 :   td->TM = shallowtrans(M);
     670        2184 :   td->PV = zero_zv(n);
     671        2184 :   gel(td->PV, 2) = Vmatrix(2, td);
     672        2184 : }
     673             : 
     674             : /* Free clones stored inside galois_test */
     675             : static void
     676        2184 : freetest(struct galois_test *td)
     677             : {
     678             :   long i;
     679       27062 :   for (i = 1; i < lg(td->PV); i++)
     680       24878 :     if (td->PV[i]) { gunclone(gel(td->PV,i)); td->PV[i] = 0; }
     681        2184 : }
     682             : 
     683             : /* Check if the integer P seen as a p-adic number is close to an integer less
     684             :  * than td->borne in absolute value */
     685             : static long
     686       64141 : padicisint(GEN P, struct galois_test *td)
     687             : {
     688       64141 :   pari_sp ltop = avma;
     689       64141 :   GEN U  = modii(P, td->ladic);
     690       64141 :   long r = cmpii(U, td->borne) <= 0 || cmpii(U, td->lborne) >= 0;
     691       64141 :   return gc_long(ltop, r);
     692             : }
     693             : 
     694             : /* Check if the permutation pf is valid according to td.
     695             :  * If not, update td to make subsequent test faster (hopefully) */
     696             : static long
     697      114058 : galois_test_perm(struct galois_test *td, GEN pf)
     698             : {
     699      114058 :   pari_sp av = avma;
     700      114058 :   long i, j, n = lg(td->L)-1;
     701      114058 :   GEN V, P = NULL;
     702      179382 :   for (i = 1; i < n; i++)
     703             :   {
     704      175749 :     long ord = td->order[i];
     705      175749 :     GEN PW = gel(td->PV, ord);
     706      175749 :     if (PW)
     707             :     {
     708      111608 :       ulong head = umael(PW,1,pf[1]);
     709     7120176 :       for (j = 2; j <= n; j++) head += umael(PW,j,pf[j]);
     710      111608 :       if (!headlongisint(head,n)) break;
     711             :     } else
     712             :     {
     713       64141 :       if (!P) P = vecpermute(td->L, pf);
     714       64141 :       V = FpV_dotproduct(gel(td->TM,ord), P, td->ladic);
     715       64141 :       if (!padicisint(V, td)) {
     716         224 :         gel(td->PV, ord) = Vmatrix(ord, td);
     717         224 :         if (DEBUGLEVEL >= 4) err_printf("M");
     718         224 :         break;
     719             :       }
     720             :     }
     721             :   }
     722      114058 :   if (i == n) return gc_long(av,1);
     723      110425 :   if (DEBUGLEVEL >= 4) err_printf("%d.", i);
     724      110425 :   if (i > 1)
     725             :   {
     726         742 :     long z = td->order[i];
     727        1526 :     for (j = i; j > 1; j--) td->order[j] = td->order[j-1];
     728         742 :     td->order[1] = z;
     729         742 :     if (DEBUGLEVEL >= 8) err_printf("%Ps", td->order);
     730             :   }
     731      110425 :   return gc_long(av,0);
     732             : }
     733             : /*Compute a*b/c when a*b will overflow*/
     734             : static long
     735           0 : muldiv(long a,long b,long c)
     736             : {
     737           0 :   return (long)((double)a*(double)b/c);
     738             : }
     739             : 
     740             : /* F = cycle decomposition of sigma,
     741             :  * B = cycle decomposition of cl(tau).
     742             :  * Check all permutations pf who can possibly correspond to tau, such that
     743             :  * tau*sigma*tau^-1 = sigma^s and tau^d = sigma^t, where d = ord cl(tau)
     744             :  * x: vector of choices,
     745             :  * G: vector allowing linear access to elts of F.
     746             :  * Choices multiple of e are not changed. */
     747             : static GEN
     748        4333 : testpermutation(GEN F, GEN B, GEN x, long s, long e, long cut,
     749             :                 struct galois_test *td)
     750             : {
     751        4333 :   pari_sp av, avm = avma;
     752             :   long a, b, c, d, n, p1, p2, p3, p4, p5, p6, l1, l2, N1, N2, R1;
     753        4333 :   long i, j, cx, hop = 0, start = 0;
     754             :   GEN pf, ar, G, W, NN, NQ;
     755             :   pari_timer ti;
     756        4333 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
     757        4333 :   a = lg(F)-1; b = lg(gel(F,1))-1;
     758        4333 :   c = lg(B)-1; d = lg(gel(B,1))-1;
     759        4333 :   n = a*b;
     760        4333 :   s = (b+s) % b;
     761        4333 :   pf = cgetg(n+1, t_VECSMALL);
     762        4333 :   av = avma;
     763        4333 :   ar = cgetg(a+2, t_VECSMALL); ar[a+1]=0;
     764        4333 :   G  = cgetg(a+1, t_VECSMALL);
     765        4333 :   W  = gel(td->PV, td->order[n]);
     766       41055 :   for (cx=1, i=1, j=1; cx <= a; cx++, i++)
     767             :   {
     768       36722 :     gel(G,cx) = gel(F, coeff(B,i,j));
     769       36722 :     if (i == d) { i = 0; j++; }
     770             :   }
     771        4333 :   NN = divis(powuu(b, c * (d - d/e)), cut);
     772        4333 :   if (DEBUGLEVEL>=4) err_printf("GaloisConj: I will try %Ps permutations\n", NN);
     773        4333 :   N1 = 1000000;
     774        4333 :   NQ = divis_rem(NN,N1,&R1);
     775        4333 :   if (abscmpiu(NQ,100000000)>0)
     776             :   {
     777           0 :     set_avma(avm);
     778           0 :     pari_warn(warner,"Combinatorics too hard: would need %Ps tests!\n"
     779             :                      "I'll skip it but you will get a partial result...",NN);
     780           0 :     return identity_perm(n);
     781             :   }
     782        4333 :   N2 = itos(NQ);
     783        5159 :   for (l2 = 0; l2 <= N2; l2++)
     784             :   {
     785        4333 :     long nbiter = (l2<N2) ? N1: R1;
     786        4333 :     if (DEBUGLEVEL >= 2 && N2) err_printf("%d%% ", muldiv(l2,100,N2));
     787    10129693 :     for (l1 = 0; l1 < nbiter; l1++)
     788             :     {
     789    10128867 :       if (start)
     790             :       {
     791    18125471 :         for (i=1, j=e; i < a;)
     792             :         {
     793    18125471 :           if ((++(x[i])) != b) break;
     794     8000937 :           x[i++] = 0;
     795     8000937 :           if (i == j) { i++; j += e; }
     796             :         }
     797             :       }
     798        4333 :       else { start=1; i = a-1; }
     799             :       /* intheadlong test: overflow in + is OK, we compute mod 2^BIL */
     800    46022151 :       for (p1 = i+1, p5 = p1%d - 1 ; p1 >= 1; p1--, p5--) /* p5 = (p1%d) - 1 */
     801             :       {
     802             :         GEN G1, G6;
     803    35893284 :         ulong V = 0;
     804    35893284 :         if (p5 == - 1) { p5 = d - 1; p6 = p1 + 1 - d; } else p6 = p1 + 1;
     805    35893284 :         G1 = gel(G,p1); G6 = gel(G,p6);
     806    35893284 :         p4 = p5 ? x[p1-1] : 0;
     807   109232354 :         for (p2 = 1+p4, p3 = 1 + x[p1]; p2 <= b; p2++)
     808             :         {
     809    73339070 :           V += umael(W,uel(G6,p3),uel(G1,p2));
     810    73339070 :           p3 += s; if (p3 > b) p3 -= b;
     811             :         }
     812    35893284 :         p3 = 1 + x[p1] - s; if (p3 <= 0) p3 += b;
     813    50186563 :         for (p2 = p4; p2 >= 1; p2--)
     814             :         {
     815    14293279 :           V += umael(W,uel(G6,p3),uel(G1,p2));
     816    14293279 :           p3 -= s; if (p3 <= 0) p3 += b;
     817             :         }
     818    35893284 :         uel(ar,p1) = uel(ar,p1+1) + V;
     819             :       }
     820    10128867 :       if (!headlongisint(uel(ar,1),n)) continue;
     821             : 
     822             :       /* intheadlong succeeds. Full computation */
     823     3498390 :       for (p1=1, p5=d; p1 <= a; p1++, p5++)
     824             :       {
     825     3384458 :         if (p5 == d) { p5 = 0; p4 = 0; } else p4 = x[p1-1];
     826     3384458 :         if (p5 == d-1) p6 = p1+1-d; else p6 = p1+1;
     827     9478385 :         for (p2 = 1+p4, p3 = 1 + x[p1]; p2 <= b; p2++)
     828             :         {
     829     6093927 :           pf[mael(G,p1,p2)] = mael(G,p6,p3);
     830     6093927 :           p3 += s; if (p3 > b) p3 -= b;
     831             :         }
     832     3384458 :         p3 = 1 + x[p1] - s; if (p3 <= 0) p3 += b;
     833     4404974 :         for (p2 = p4; p2 >= 1; p2--)
     834             :         {
     835     1020516 :           pf[mael(G,p1,p2)] = mael(G,p6,p3);
     836     1020516 :           p3 -= s; if (p3 <= 0) p3 += b;
     837             :         }
     838             :       }
     839      113932 :       if (galois_test_perm(td, pf))
     840             :       {
     841        3507 :         if (DEBUGLEVEL >= 1)
     842             :         {
     843           0 :           GEN nb = addis(mulss(l2,N1),l1);
     844           0 :           timer_printf(&ti, "testpermutation(%Ps)", nb);
     845           0 :           if (DEBUGLEVEL >= 2 && hop)
     846           0 :             err_printf("GaloisConj: %d hop over %Ps iterations\n", hop, nb);
     847             :         }
     848        3507 :         set_avma(av); return pf;
     849             :       }
     850      110425 :       hop++;
     851             :     }
     852             :   }
     853         826 :   if (DEBUGLEVEL >= 1)
     854             :   {
     855           0 :     timer_printf(&ti, "testpermutation(%Ps)", NN);
     856           0 :     if (DEBUGLEVEL >= 2 && hop)
     857           0 :       err_printf("GaloisConj: %d hop over %Ps iterations\n", hop, NN);
     858             :   }
     859         826 :   return gc_NULL(avm);
     860             : }
     861             : 
     862             : /* List of subgroups of (Z/mZ)^* whose order divide o, and return the list
     863             :  * of their elements, sorted by increasing order */
     864             : static GEN
     865         525 : listznstarelts(long m, long o)
     866             : {
     867         525 :   pari_sp av = avma;
     868             :   GEN L, zn, zns;
     869             :   long i, phi, ind, l;
     870         525 :   if (m == 2) retmkvec(mkvecsmall(1));
     871         511 :   zn = znstar(stoi(m));
     872         511 :   phi = itos(gel(zn,1));
     873         511 :   o = ugcd(o, phi); /* do we impose this on input ? */
     874         511 :   zns = znstar_small(zn);
     875         511 :   L = cgetg(o+1, t_VEC);
     876        1547 :   for (i=1,ind = phi; ind; ind -= phi/o, i++) /* by *decreasing* exact index */
     877        1036 :     gel(L,i) = subgrouplist(gel(zn,2), mkvec(utoipos(ind)));
     878         511 :   L = shallowconcat1(L); l = lg(L);
     879        1526 :   for (i = 1; i < l; i++) gel(L,i) = znstar_hnf_elts(zns, gel(L,i));
     880         511 :   return gerepilecopy(av, L);
     881             : }
     882             : 
     883             : /* A sympol is a symmetric polynomial
     884             :  *
     885             :  * Currently sympol are couple of t_VECSMALL [v,w]
     886             :  * v[1]...v[k], w[1]...w[k]  represent the polynomial sum(i=1,k,v[i]*s_w[i])
     887             :  * where s_i(X_1,...,X_n) = sum(j=1,n,X_j^i) */
     888             : 
     889             : static GEN
     890        5034 : Flm_newtonsum(GEN M, ulong e, ulong p)
     891             : {
     892        5034 :   long f = lg(M), g = lg(gel(M,1)), i, j;
     893        5034 :   GEN NS = cgetg(f, t_VECSMALL);
     894       28383 :   for(i=1; i<f; i++)
     895             :   {
     896       23349 :     ulong s = 0;
     897       23349 :     GEN Mi = gel(M,i);
     898      102223 :     for(j = 1; j < g; j++)
     899       78874 :       s = Fl_add(s, Fl_powu(uel(Mi,j), e, p), p);
     900       23349 :     uel(NS,i) = s;
     901             :   }
     902        5034 :   return NS;
     903             : }
     904             : 
     905             : static GEN
     906        3437 : Flv_sympol_eval(GEN v, GEN NS, ulong p)
     907             : {
     908        3437 :   pari_sp av = avma;
     909        3437 :   long i, l = lg(v);
     910        3437 :   GEN S = Flv_Fl_mul(gel(NS,1), uel(v,1), p);
     911        3731 :   for (i=2; i<l; i++)
     912         294 :     if (v[i]) S = Flv_add(S, Flv_Fl_mul(gel(NS,i), uel(v,i), p), p);
     913        3437 :   return gerepileuptoleaf(av, S);
     914             : }
     915             : 
     916             : static GEN
     917        3437 : sympol_eval_newtonsum(long e, GEN O, GEN mod)
     918             : {
     919        3437 :   long f = lg(O), g = lg(gel(O,1)), i, j;
     920        3437 :   GEN PL = cgetg(f, t_COL);
     921       20601 :   for(i=1; i<f; i++)
     922             :   {
     923       17164 :     pari_sp av = avma;
     924       17164 :     GEN s = gen_0;
     925       66318 :     for(j=1; j<g; j++) s = addii(s, Fp_powu(gmael(O,i,j), e, mod));
     926       17164 :     gel(PL,i) = gerepileuptoint(av, remii(s,mod));
     927             :   }
     928        3437 :   return PL;
     929             : }
     930             : 
     931             : static GEN
     932        3374 : sympol_eval(GEN sym, GEN O, GEN mod)
     933             : {
     934        3374 :   pari_sp av = avma;
     935             :   long i;
     936        3374 :   GEN v = gel(sym,1), w = gel(sym,2);
     937        3374 :   GEN S = gen_0;
     938        6979 :   for (i=1; i<lg(v); i++)
     939        3605 :     if (v[i]) S = gadd(S, gmulsg(v[i],  sympol_eval_newtonsum(w[i], O, mod)));
     940        3374 :   return gerepileupto(av, S);
     941             : }
     942             : 
     943             : /* Let sigma be an automorphism of L (as a polynomial with rational coefs)
     944             :  * Let 'sym' be a symmetric polynomial defining alpha in L.
     945             :  * We have alpha = sym(x,sigma(x),,,sigma^(g-1)(x)). Compute alpha mod p */
     946             : static GEN
     947        2163 : sympol_aut_evalmod(GEN sym, long g, GEN sigma, GEN Tp, GEN p)
     948             : {
     949        2163 :   pari_sp ltop=avma;
     950        2163 :   long i, j, npows = brent_kung_optpow(degpol(Tp)-1, g-1, 1);
     951        2163 :   GEN s, f, pows, v = zv_to_ZV(gel(sym,1)), w = zv_to_ZV(gel(sym,2));
     952        2163 :   sigma = RgX_to_FpX(sigma, p);
     953        2163 :   pows  = FpXQ_powers(sigma,npows,Tp,p);
     954        2163 :   f = pol_x(varn(sigma));
     955        2163 :   s = pol_0(varn(sigma));
     956        8911 :   for(i=1; i<=g;i++)
     957             :   {
     958        6748 :     if (i > 1) f = FpX_FpXQV_eval(f,pows,Tp,p);
     959       13874 :     for(j=1; j<lg(v); j++)
     960        7126 :       s = FpX_add(s, FpX_Fp_mul(FpXQ_pow(f,gel(w,j),Tp,p),gel(v,j),p),p);
     961             :   }
     962        2163 :   return gerepileupto(ltop, s);
     963             : }
     964             : 
     965             : /* Let Sp be as computed with sympol_aut_evalmod
     966             :  * Let Tmod be the factorisation of T mod p.
     967             :  * Return the factorisation of the minimal polynomial of S mod p w.r.t. Tmod */
     968             : static GEN
     969        2163 : fixedfieldfactmod(GEN Sp, GEN p, GEN Tmod)
     970             : {
     971        2163 :   long i, l = lg(Tmod);
     972        2163 :   GEN F = cgetg(l,t_VEC);
     973        8302 :   for(i=1; i<l; i++)
     974             :   {
     975        6139 :     GEN Ti = gel(Tmod,i);
     976        6139 :     gel(F,i) = FpXQ_minpoly(FpX_rem(Sp,Ti,p), Ti,p);
     977             :   }
     978        2163 :   return F;
     979             : }
     980             : 
     981             : static GEN
     982        3374 : fixedfieldsurmer(ulong l, GEN NS, GEN W)
     983             : {
     984        3374 :   const long step=3;
     985        3374 :   long i, j, n = lg(W)-1, m = 1L<<((n-1)<<1);
     986        3374 :   GEN sym = cgetg(n+1,t_VECSMALL);
     987        3605 :   for (j=1;j<n;j++) sym[j] = step;
     988        3374 :   sym[n] = 0;
     989        3374 :   if (DEBUGLEVEL>=4) err_printf("FixedField: Weight: %Ps\n",W);
     990        3437 :   for (i=0; i<m; i++)
     991             :   {
     992        3437 :     pari_sp av = avma;
     993             :     GEN L;
     994        3668 :     for (j=1; sym[j]==step; j++) sym[j]=0;
     995        3437 :     sym[j]++;
     996        3437 :     if (DEBUGLEVEL>=6) err_printf("FixedField: Sym: %Ps\n",sym);
     997        3437 :     L = Flv_sympol_eval(sym, NS, l);
     998        3437 :     if (!vecsmall_is1to1(L)) { set_avma(av); continue; }
     999        3374 :     return mkvec2(sym,W);
    1000             :   }
    1001           0 :   return NULL;
    1002             : }
    1003             : 
    1004             : /*Check whether the line of NS are pair-wise distinct.*/
    1005             : static long
    1006        3605 : sympol_is1to1_lg(GEN NS, long n)
    1007             : {
    1008        3605 :   long i, j, k, l = lgcols(NS);
    1009       19502 :   for (i=1; i<l; i++)
    1010       92253 :     for(j=i+1; j<l; j++)
    1011             :     {
    1012       78743 :       for(k=1; k<n; k++)
    1013       78512 :         if (mael(NS,k,j)!=mael(NS,k,i)) break;
    1014       76356 :       if (k>=n) return 0;
    1015             :     }
    1016        3374 :   return 1;
    1017             : }
    1018             : 
    1019             : /* Let O a set of orbits of roots (see fixedfieldorbits) modulo mod,
    1020             :  * l | mod and p two prime numbers. Return a vector [sym,s,P] where:
    1021             :  * sym is a sympol, s is the set of images of sym on O and
    1022             :  * P is the polynomial with roots s. */
    1023             : static GEN
    1024        3374 : fixedfieldsympol(GEN O, ulong l)
    1025             : {
    1026        3374 :   pari_sp ltop=avma;
    1027        3374 :   const long n=(BITS_IN_LONG>>1)-1;
    1028        3374 :   GEN NS = cgetg(n+1,t_MAT), sym = NULL, W = cgetg(n+1,t_VECSMALL);
    1029        3374 :   long i, e=1;
    1030        3374 :   if (DEBUGLEVEL>=4)
    1031           0 :     err_printf("FixedField: Size: %ldx%ld\n",lg(O)-1,lg(gel(O,1))-1);
    1032        3374 :   O = ZM_to_Flm(O,l);
    1033        6979 :   for (i=1; !sym && i<=n; i++)
    1034             :   {
    1035        3605 :     GEN L = Flm_newtonsum(O, e++, l);
    1036        3605 :     if (lg(O)>2)
    1037        4971 :       while (vecsmall_isconst(L)) L = Flm_newtonsum(O, e++, l);
    1038        3605 :     W[i] = e-1; gel(NS,i) = L;
    1039        3605 :     if (sympol_is1to1_lg(NS,i+1))
    1040        3374 :       sym = fixedfieldsurmer(l,NS,vecsmall_shorten(W,i));
    1041             :   }
    1042        3374 :   if (!sym) pari_err_BUG("fixedfieldsympol [p too small]");
    1043        3374 :   if (DEBUGLEVEL>=2) err_printf("FixedField: Found: %Ps\n",gel(sym,1));
    1044        3374 :   return gerepilecopy(ltop,sym);
    1045             : }
    1046             : 
    1047             : /* Let O a set of orbits as indices and L the corresponding roots.
    1048             :  * Return the set of orbits as roots. */
    1049             : static GEN
    1050        3374 : fixedfieldorbits(GEN O, GEN L)
    1051             : {
    1052        3374 :   GEN S = cgetg(lg(O), t_MAT);
    1053             :   long i;
    1054       19222 :   for (i = 1; i < lg(O); i++) gel(S,i) = vecpermute(L, gel(O,i));
    1055        3374 :   return S;
    1056             : }
    1057             : 
    1058             : static GEN
    1059         875 : fixedfieldinclusion(GEN O, GEN PL)
    1060             : {
    1061         875 :   long i, j, f = lg(O)-1, g = lg(gel(O,1))-1;
    1062         875 :   GEN S = cgetg(f*g + 1, t_COL);
    1063        6167 :   for (i = 1; i <= f; i++)
    1064             :   {
    1065        5292 :     GEN Oi = gel(O,i);
    1066       22834 :     for (j = 1; j <= g; j++) gel(S, Oi[j]) = gel(PL, i);
    1067             :   }
    1068         875 :   return S;
    1069             : }
    1070             : 
    1071             : /* Polynomial attached to a vector of conjugates. Not stack clean */
    1072             : static GEN
    1073       21035 : vectopol(GEN v, GEN M, GEN den , GEN mod, GEN mod2, long x)
    1074             : {
    1075       21035 :   long l = lg(v)+1, i;
    1076       21035 :   GEN z = cgetg(l,t_POL);
    1077       21035 :   z[1] = evalsigne(1)|evalvarn(x);
    1078      257852 :   for (i=2; i<l; i++)
    1079      236817 :     gel(z,i) = gdiv(centermodii(ZMrow_ZC_mul(M,v,i-1), mod, mod2), den);
    1080       21035 :   return normalizepol_lg(z, l);
    1081             : }
    1082             : 
    1083             : /* Polynomial associate to a permutation of the roots. Not stack clean */
    1084             : static GEN
    1085       19516 : permtopol(GEN p, GEN L, GEN M, GEN den, GEN mod, GEN mod2, long x)
    1086             : {
    1087       19516 :   if (lg(p) != lg(L)) pari_err_TYPE("permtopol [permutation]", p);
    1088       19516 :   return vectopol(vecpermute(L,p), M, den, mod, mod2, x);
    1089             : }
    1090             : 
    1091             : static GEN
    1092        1183 : galoisvecpermtopol(GEN gal, GEN vec, GEN mod, GEN mod2)
    1093             : {
    1094        1183 :   long i, l = lg(vec);
    1095        1183 :   long v = varn(gal_get_pol(gal));
    1096        1183 :   GEN L = gal_get_roots(gal);
    1097        1183 :   GEN M = gal_get_invvdm(gal);
    1098        1183 :   GEN P = cgetg(l, t_MAT);
    1099        6797 :   for (i=1; i<l; i++)
    1100        5614 :     gel(P, i) = vecpermute(L,gel(vec,i));
    1101        1183 :   P = RgM_to_RgXV(FpM_center(FpM_mul(M, P, mod), mod, mod2), v);
    1102        1183 :   return gdiv(P, gal_get_den(gal));
    1103             : }
    1104             : 
    1105             : static void
    1106        1841 : notgalois(long p, struct galois_analysis *ga)
    1107             : {
    1108        1841 :   if (DEBUGLEVEL >= 2) err_printf("GaloisAnalysis:non Galois for p=%ld\n", p);
    1109        1841 :   ga->p = p;
    1110        1841 :   ga->deg = 0;
    1111        1841 : }
    1112             : 
    1113             : /*Gather information about the group*/
    1114             : static long
    1115        4942 : init_group(long n, long np, GEN Fp, GEN Fe, long *porder)
    1116             : {
    1117        4942 :   const long prim_nonwss_orders[] = { 36,48,56,60,75,80,196,200 };
    1118        4942 :   long i, phi_order = 1, order = 1, group = 0;
    1119             :   ulong p;
    1120             : 
    1121             :  /* non-WSS groups of this order? */
    1122       44331 :   for (i=0; i < (long)numberof(prim_nonwss_orders); i++)
    1123       39410 :     if (n % prim_nonwss_orders[i] == 0) { group |= ga_non_wss; break; }
    1124        4942 :   if (np == 2 && Fp[2] == 3 && Fe[2] == 1 && Fe[1] > 2) group |= ga_ext_2;
    1125             : 
    1126        7350 :   for (i = np; i > 0; i--)
    1127             :   {
    1128        5908 :     long p = Fp[i];
    1129        5908 :     if (phi_order % p == 0) { group |= ga_all_normal; break; }
    1130        4984 :     order *= p; phi_order *= p-1;
    1131        4984 :     if (Fe[i] > 1) break;
    1132             :   }
    1133        4942 :   if (uisprimepower(n, &p) || n == 135) group |= ga_all_nilpotent;
    1134        4942 :   if (n <= 104) group |= ga_easy; /* no need to use polynomial algo */
    1135        4942 :   *porder = order; return group;
    1136             : }
    1137             : 
    1138             : /*is a "better" than b ? (if so, update karma) */
    1139             : static int
    1140       22449 : improves(long a, long b, long plift, long p, long n, long *karma)
    1141             : {
    1142       22449 :   if (!plift || a > b) { *karma = ugcd(p-1,n); return 1; }
    1143       21182 :   if (a == b) {
    1144       18753 :     long k = ugcd(p-1,n);
    1145       18753 :     if (k > *karma) { *karma = k; return 1; }
    1146             :   }
    1147       19208 :   return 0; /* worse */
    1148             : }
    1149             : 
    1150             : /* return 0 if not galois or not wss */
    1151             : static int
    1152        4942 : galoisanalysis(GEN T, struct galois_analysis *ga, long calcul_l, GEN bad)
    1153             : {
    1154        4942 :   pari_sp ltop = avma, av;
    1155        4942 :   long group, linf, n, p, i, karma = 0;
    1156             :   GEN F, Fp, Fe, Fpe, O;
    1157             :   long np, order, plift, nbmax, nbtest, deg;
    1158             :   pari_timer ti;
    1159             :   forprime_t S;
    1160        4942 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
    1161        4942 :   n = degpol(T);
    1162        4942 :   O = zero_zv(n);
    1163        4942 :   F = factoru_pow(n);
    1164        4942 :   Fp = gel(F,1); np = lg(Fp)-1;
    1165        4942 :   Fe = gel(F,2);
    1166        4942 :   Fpe= gel(F,3);
    1167        4942 :   group = init_group(n, np, Fp, Fe, &order);
    1168             : 
    1169             :   /*Now we study the orders of the Frobenius elements*/
    1170        4942 :   deg = Fp[np]; /* largest prime | n */
    1171        4942 :   plift = 0;
    1172        4942 :   nbtest = 0;
    1173        4942 :   nbmax = 8+(n>>1);
    1174        4942 :   u_forprime_init(&S, n*maxss(expu(n)-3, 2), ULONG_MAX);
    1175        4942 :   av = avma;
    1176       44359 :   while (!plift || (nbtest < nbmax && (nbtest <=8 || order < (n>>1)))
    1177        3122 :                 || (n == 24 && O[6] == 0 && O[4] == 0)
    1178        3122 :                 || ((group&ga_non_wss) && order == Fp[np]))
    1179             :   {
    1180       41237 :     long d, o, norm_o = 1;
    1181             :     GEN D, Tp;
    1182             : 
    1183       41237 :     if ((group&ga_non_wss) && nbtest >= 3*nbmax) break; /* in all cases */
    1184       41237 :     nbtest++; set_avma(av);
    1185       41237 :     p = u_forprime_next(&S);
    1186       41237 :     if (!p) pari_err_OVERFLOW("galoisanalysis [ran out of primes]");
    1187       43925 :     if (bad && dvdiu(bad, p)) continue;
    1188       41237 :     Tp = ZX_to_Flx(T,p);
    1189       41237 :     if (!Flx_is_squarefree(Tp,p)) { if (!--nbtest) nbtest = 1; continue; }
    1190             : 
    1191       38549 :     D = Flx_nbfact_by_degree(Tp, &d, p);
    1192       38549 :     o = n / d; /* d factors, all should have degree o */
    1193       38549 :     if (D[o] != d) { notgalois(p, ga); return gc_bool(ltop,0); }
    1194             : 
    1195       36729 :     if (!O[o]) O[o] = p;
    1196       36729 :     if (o % deg) goto ga_end; /* NB: deg > 1 */
    1197       26278 :     if ((group&ga_all_normal) && o < order) goto ga_end;
    1198             : 
    1199             :     /*Frob_p has order o > 1, find a power which generates a normal subgroup*/
    1200       26124 :     if (o * Fp[1] >= n)
    1201       15939 :       norm_o = o; /*subgroups of smallest index are normal*/
    1202             :     else
    1203             :     {
    1204       11361 :       for (i = np; i > 0; i--)
    1205             :       {
    1206       11361 :         if (o % Fpe[i]) break;
    1207        1176 :         norm_o *= Fpe[i];
    1208             :       }
    1209             :     }
    1210             :     /* Frob_p^(o/norm_o) generates a normal subgroup of order norm_o */
    1211       26124 :     if (norm_o != 1)
    1212             :     {
    1213       17115 :       if (!(group&ga_all_normal) || o > order)
    1214        3675 :         karma = ugcd(p-1,n);
    1215       13440 :       else if (!improves(norm_o, deg, plift,p,n, &karma)) goto ga_end;
    1216             :       /* karma0=0, deg0<=norm_o -> the first improves() returns 1 */
    1217        5761 :       deg = norm_o; group |= ga_all_normal; /* STORE */
    1218             :     }
    1219        9009 :     else if (group&ga_all_normal) goto ga_end;
    1220        9009 :     else if (!improves(o, order, plift,p,n, &karma)) goto ga_end;
    1221             : 
    1222        6916 :     order = o; plift = p; /* STORE */
    1223       36729 :     ga_end:
    1224       36729 :     if (DEBUGLEVEL >= 5)
    1225           0 :       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);
    1226             :   }
    1227             :   /* To avoid looping on non-wss group.
    1228             :    * TODO: check for large groups. Would it be better to disable this check if
    1229             :    * we are in a good case (ga_all_normal && !(ga_ext_2) (e.g. 60)) ?*/
    1230        3122 :   ga->p = plift;
    1231        3122 :   if (!plift || ((group&ga_non_wss) && order == Fp[np]))
    1232             :   {
    1233           0 :     pari_warn(warner,"Galois group almost certainly not weakly super solvable");
    1234           0 :     return 0;
    1235             :   }
    1236        3122 :   linf = 2*n*usqrt(n);
    1237        3122 :   if (calcul_l && O[1] <= linf)
    1238             :   {
    1239             :     pari_sp av2;
    1240             :     forprime_t S2;
    1241             :     ulong p;
    1242        1449 :     u_forprime_init(&S2, linf+1,ULONG_MAX);
    1243        1449 :     av2 = avma;
    1244       39095 :     while ((p = u_forprime_next(&S2)))
    1245             :     { /*find a totally split prime l > linf*/
    1246       39095 :       GEN Tp = ZX_to_Flx(T, p);
    1247       39095 :       long nb = Flx_nbroots(Tp, p);
    1248       39095 :       if (nb == n) { O[1] = p; break; }
    1249       37667 :       if (nb && Flx_is_squarefree(Tp,p)) { notgalois(p,ga); return gc_bool(ltop,0); }
    1250       37646 :       set_avma(av2);
    1251             :     }
    1252        1428 :     if (!p) pari_err_OVERFLOW("galoisanalysis [ran out of primes]");
    1253             :   }
    1254        3101 :   ga->group = group;
    1255        3101 :   ga->deg = deg;
    1256        3101 :   ga->ord = order;
    1257        3101 :   ga->l  = O[1];
    1258        3101 :   ga->p4 = n >= 4 ? O[4] : 0;
    1259        3101 :   if (DEBUGLEVEL >= 4)
    1260           0 :     err_printf("GaloisAnalysis:p=%ld l=%ld group=%ld deg=%ld ord=%ld\n",
    1261           0 :                plift, O[1], group, deg, order);
    1262        3101 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "galoisanalysis()");
    1263        3101 :   return gc_bool(ltop,1);
    1264             : }
    1265             : 
    1266             : static GEN
    1267          42 : a4galoisgen(struct galois_test *td)
    1268             : {
    1269          42 :   const long n = 12;
    1270          42 :   pari_sp ltop = avma, av, av2;
    1271          42 :   long i, j, k, N, hop = 0;
    1272             :   GEN MT, O,O1,O2,O3, ar, mt, t, u, res, orb, pft, pfu, pfv;
    1273             : 
    1274          42 :   res = cgetg(3, t_VEC);
    1275          42 :   pft = cgetg(n+1, t_VECSMALL);
    1276          42 :   pfu = cgetg(n+1, t_VECSMALL);
    1277          42 :   pfv = cgetg(n+1, t_VECSMALL);
    1278          42 :   gel(res,1) = mkvec3(pft,pfu,pfv);
    1279          42 :   gel(res,2) = mkvecsmall3(2,2,3);
    1280          42 :   av = avma;
    1281          42 :   ar = cgetg(5, t_VECSMALL);
    1282          42 :   mt = gel(td->PV, td->order[n]);
    1283          42 :   t = identity_perm(n) + 1; /* Sorry for this hack */
    1284          42 :   u = cgetg(n+1, t_VECSMALL) + 1; /* too lazy to correct */
    1285          42 :   MT = cgetg(n+1, t_MAT);
    1286         546 :   for (j = 1; j <= n; j++) gel(MT,j) = cgetg(n+1, t_VECSMALL);
    1287         546 :   for (j = 1; j <= n; j++)
    1288        3276 :     for (i = 1; i < j; i++)
    1289        2772 :       ucoeff(MT,i,j) = ucoeff(MT,j,i) = ucoeff(mt,i,j)+ucoeff(mt,j,i);
    1290             :   /* MT(i,i) unused */
    1291             : 
    1292          42 :   av2 = avma;
    1293             :   /* N = itos(gdiv(mpfact(n), mpfact(n >> 1))) >> (n >> 1); */
    1294             :   /* n = 2k = 12; N = (2k)! / (k! * 2^k) = 10395 */
    1295          42 :   N = 10395;
    1296          42 :   if (DEBUGLEVEL>=4) err_printf("A4GaloisConj: will test %ld permutations\n", N);
    1297          42 :   uel(ar,4) = umael(MT,11,12);
    1298          42 :   uel(ar,3) = uel(ar,4) + umael(MT,9,10);
    1299          42 :   uel(ar,2) = uel(ar,3) + umael(MT,7,8);
    1300          42 :   uel(ar,1) = uel(ar,2) + umael(MT,5,6);
    1301      102569 :   for (i = 0; i < N; i++)
    1302             :   {
    1303             :     long g;
    1304      102569 :     if (i)
    1305             :     {
    1306      102527 :       long a, x = i, y = 1;
    1307      144536 :       do { y += 2; a = x%y; x = x/y; } while (!a);
    1308      102527 :       switch (y)
    1309             :       {
    1310       68368 :       case 3:
    1311       68368 :         lswap(t[2], t[2-a]);
    1312       68368 :         break;
    1313       27341 :       case 5:
    1314       27341 :         x = t[0]; t[0] = t[2]; t[2] = t[1]; t[1] = x;
    1315       27341 :         lswap(t[4], t[4-a]);
    1316       27341 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1317       27341 :         break;
    1318        5865 :       case 7:
    1319        5865 :         x = t[0]; t[0] = t[4]; t[4] = t[3]; t[3] = t[1]; t[1] = t[2]; t[2] = x;
    1320        5865 :         lswap(t[6], t[6-a]);
    1321        5865 :         uel(ar,2) = uel(ar,3) + umael(MT,t[6],t[7]);
    1322        5865 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1323        5865 :         break;
    1324         874 :       case 9:
    1325         874 :         x = t[0]; t[0] = t[6]; t[6] = t[5]; t[5] = t[3]; t[3] = x;
    1326         874 :         lswap(t[1], t[4]);
    1327         874 :         lswap(t[8], t[8-a]);
    1328         874 :         uel(ar,3) = uel(ar,4) + umael(MT,t[8],t[9]);
    1329         874 :         uel(ar,2) = uel(ar,3) + umael(MT,t[6],t[7]);
    1330         874 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1331         874 :         break;
    1332          79 :       case 11:
    1333          79 :         x = t[0]; t[0] = t[8]; t[8] = t[7]; t[7] = t[5]; t[5] = t[1];
    1334          79 :         t[1] = t[6]; t[6] = t[3]; t[3] = t[2]; t[2] = t[4]; t[4] = x;
    1335          79 :         lswap(t[10], t[10-a]);
    1336          79 :         uel(ar,4) = umael(MT,t[10],t[11]);
    1337          79 :         uel(ar,3) = uel(ar,4) + umael(MT,t[8],t[9]);
    1338          79 :         uel(ar,2) = uel(ar,3) + umael(MT,t[6],t[7]);
    1339          79 :         uel(ar,1) = uel(ar,2) + umael(MT,t[4],t[5]);
    1340             :       }
    1341          42 :     }
    1342      102569 :     g = uel(ar,1)+umael(MT,t[0],t[1])+umael(MT,t[2],t[3]);
    1343      102569 :     if (headlongisint(g,n))
    1344             :     {
    1345         294 :       for (k = 0; k < n; k += 2)
    1346             :       {
    1347         252 :         pft[t[k]] = t[k+1];
    1348         252 :         pft[t[k+1]] = t[k];
    1349             :       }
    1350          42 :       if (galois_test_perm(td, pft)) break;
    1351           0 :       hop++;
    1352             :     }
    1353      102527 :     set_avma(av2);
    1354             :   }
    1355          42 :   if (DEBUGLEVEL >= 1 && hop)
    1356           0 :     err_printf("A4GaloisConj: %ld hop over %ld iterations\n", hop, N);
    1357          42 :   if (i == N) return gc_NULL(ltop);
    1358             :   /* N = itos(gdiv(mpfact(n >> 1), mpfact(n >> 2))) >> 1; */
    1359          42 :   N = 60;
    1360          42 :   if (DEBUGLEVEL >= 4) err_printf("A4GaloisConj: sigma=%Ps \n", pft);
    1361         168 :   for (k = 0; k < n; k += 4)
    1362             :   {
    1363         126 :     u[k+3] = t[k+3];
    1364         126 :     u[k+2] = t[k+1];
    1365         126 :     u[k+1] = t[k+2];
    1366         126 :     u[k]   = t[k];
    1367             :   }
    1368        2018 :   for (i = 0; i < N; i++)
    1369             :   {
    1370        2018 :     ulong g = 0;
    1371        2018 :     if (i)
    1372             :     {
    1373        1976 :       long a, x = i, y = -2;
    1374        3110 :       do { y += 4; a = x%y; x = x/y; } while (!a);
    1375        1976 :       lswap(u[0],u[2]);
    1376        1976 :       switch (y)
    1377             :       {
    1378         988 :       case 2:
    1379         988 :         break;
    1380         842 :       case 6:
    1381         842 :         lswap(u[4],u[6]);
    1382         842 :         if (!(a & 1))
    1383             :         {
    1384         341 :           a = 4 - (a>>1);
    1385         341 :           lswap(u[6], u[a]);
    1386         341 :           lswap(u[4], u[a-2]);
    1387             :         }
    1388         842 :         break;
    1389         146 :       case 10:
    1390         146 :         x = u[6];
    1391         146 :         u[6] = u[3];
    1392         146 :         u[3] = u[2];
    1393         146 :         u[2] = u[4];
    1394         146 :         u[4] = u[1];
    1395         146 :         u[1] = u[0];
    1396         146 :         u[0] = x;
    1397         146 :         if (a >= 3) a += 2;
    1398         146 :         a = 8 - a;
    1399         146 :         lswap(u[10],u[a]);
    1400         146 :         lswap(u[8], u[a-2]);
    1401         146 :         break;
    1402             :       }
    1403          42 :     }
    1404       14126 :     for (k = 0; k < n; k += 2) g += mael(MT,u[k],u[k+1]);
    1405        2018 :     if (headlongisint(g,n))
    1406             :     {
    1407         294 :       for (k = 0; k < n; k += 2)
    1408             :       {
    1409         252 :         pfu[u[k]] = u[k+1];
    1410         252 :         pfu[u[k+1]] = u[k];
    1411             :       }
    1412          42 :       if (galois_test_perm(td, pfu)) break;
    1413           0 :       hop++;
    1414             :     }
    1415        1976 :     set_avma(av2);
    1416             :   }
    1417          42 :   if (i == N) return gc_NULL(ltop);
    1418          42 :   if (DEBUGLEVEL >= 1 && hop)
    1419           0 :     err_printf("A4GaloisConj: %ld hop over %ld iterations\n", hop, N);
    1420          42 :   if (DEBUGLEVEL >= 4) err_printf("A4GaloisConj: tau=%Ps \n", pfu);
    1421          42 :   set_avma(av2);
    1422          42 :   orb = mkvec2(pft,pfu);
    1423          42 :   O = vecperm_orbits(orb, 12);
    1424          42 :   if (DEBUGLEVEL >= 4) {
    1425           0 :     err_printf("A4GaloisConj: orb=%Ps\n", orb);
    1426           0 :     err_printf("A4GaloisConj: O=%Ps \n", O);
    1427             :   }
    1428          42 :   av2 = avma;
    1429          42 :   O1 = gel(O,1); O2 = gel(O,2); O3 = gel(O,3);
    1430          63 :   for (j = 0; j < 2; j++)
    1431             :   {
    1432          63 :     pfv[O1[1]] = O2[1];
    1433          63 :     pfv[O1[2]] = O2[3+j];
    1434          63 :     pfv[O1[3]] = O2[4 - (j << 1)];
    1435          63 :     pfv[O1[4]] = O2[2+j];
    1436         203 :     for (i = 0; i < 4; i++)
    1437             :     {
    1438         182 :       ulong g = 0;
    1439         182 :       switch (i)
    1440             :       {
    1441          63 :       case 0: break;
    1442          56 :       case 1: lswap(O3[1], O3[2]); lswap(O3[3], O3[4]); break;
    1443          42 :       case 2: lswap(O3[1], O3[4]); lswap(O3[2], O3[3]); break;
    1444          21 :       case 3: lswap(O3[1], O3[2]); lswap(O3[3], O3[4]); break;
    1445             :       }
    1446         182 :       pfv[O2[1]]          = O3[1];
    1447         182 :       pfv[O2[3+j]]        = O3[4-j];
    1448         182 :       pfv[O2[4 - (j<<1)]] = O3[2 + (j<<1)];
    1449         182 :       pfv[O2[2+j]]        = O3[3-j];
    1450         182 :       pfv[O3[1]]          = O1[1];
    1451         182 :       pfv[O3[4-j]]        = O1[2];
    1452         182 :       pfv[O3[2 + (j<<1)]] = O1[3];
    1453         182 :       pfv[O3[3-j]]        = O1[4];
    1454        2366 :       for (k = 1; k <= n; k++) g += mael(mt,k,pfv[k]);
    1455         182 :       if (headlongisint(g,n) && galois_test_perm(td, pfv))
    1456             :       {
    1457          42 :         set_avma(av);
    1458          42 :         if (DEBUGLEVEL >= 1)
    1459           0 :           err_printf("A4GaloisConj: %ld hop over %d iterations max\n",
    1460             :                      hop, 10395 + 68);
    1461          42 :         return res;
    1462             :       }
    1463         140 :       hop++; set_avma(av2);
    1464             :     }
    1465             :   }
    1466           0 :   return gc_NULL(ltop);
    1467             : }
    1468             : 
    1469             : /* S4 */
    1470             : static void
    1471         217 : s4makelift(GEN u, struct galois_lift *gl, GEN liftpow)
    1472             : {
    1473         217 :   GEN s = automorphismlift(u, gl);
    1474             :   long i;
    1475         217 :   gel(liftpow,1) = s;
    1476        4991 :   for (i = 2; i < lg(liftpow); i++)
    1477        4774 :     gel(liftpow,i) = FpXQ_mul(gel(liftpow,i-1), s, gl->TQ, gl->Q);
    1478         217 : }
    1479             : static long
    1480        3481 : s4test(GEN u, GEN liftpow, struct galois_lift *gl, GEN phi)
    1481             : {
    1482        3481 :   pari_sp av = avma;
    1483             :   GEN res, Q, Q2;
    1484        3481 :   long bl, i, d = lg(u)-2;
    1485             :   pari_timer ti;
    1486        3481 :   if (DEBUGLEVEL >= 6) timer_start(&ti);
    1487        3481 :   if (!d) return 0;
    1488        3481 :   Q = gl->Q; Q2 = shifti(Q,-1);
    1489        3481 :   res = gel(u,2);
    1490       81828 :   for (i = 1; i < d; i++)
    1491       78347 :     if (lg(gel(liftpow,i))>2)
    1492       78347 :       res = addii(res, mulii(gmael(liftpow,i,2), gel(u,i+2)));
    1493        3481 :   res = remii(res,Q);
    1494        3481 :   if (gl->den != gen_1) res = mulii(res, gl->den);
    1495        3481 :   res = centermodii(res, Q,Q2);
    1496        3481 :   if (abscmpii(res, gl->gb->bornesol) > 0) return gc_long(av,0);
    1497         393 :   res = scalar_ZX_shallow(gel(u,2),varn(u));
    1498        8948 :   for (i = 1; i < d ; i++)
    1499        8555 :     if (lg(gel(liftpow,i))>2)
    1500        8555 :       res = ZX_add(res, ZX_Z_mul(gel(liftpow,i), gel(u,i+2)));
    1501         393 :   res = FpX_red(res, Q);
    1502         393 :   if (gl->den != gen_1) res = FpX_Fp_mul(res, gl->den, Q);
    1503         393 :   res = FpX_center_i(res, Q, shifti(Q,-1));
    1504         393 :   bl = poltopermtest(res, gl, phi);
    1505         393 :   if (DEBUGLEVEL >= 6) timer_printf(&ti, "s4test()");
    1506         393 :   return gc_long(av,bl);
    1507             : }
    1508             : 
    1509             : static GEN
    1510         651 : aux(long a, long b, GEN T, GEN M, GEN p, GEN *pu)
    1511             : {
    1512         651 :   *pu = FpX_mul(gel(T,b), gel(T,a),p);
    1513        2604 :   return FpX_chinese_coprime(gmael(M,a,b), gmael(M,b,a),
    1514         651 :                              gel(T,b), gel(T,a), *pu, p);
    1515             : }
    1516             : 
    1517             : static GEN
    1518         217 : s4releveauto(GEN misom,GEN Tmod,GEN Tp,GEN p,long a1,long a2,long a3,long a4,long a5,long a6)
    1519             : {
    1520         217 :   pari_sp av = avma;
    1521             :   GEN u4,u5;
    1522             :   GEN pu1, pu2, pu3, pu4;
    1523         217 :   GEN u1 = aux(a1, a2, Tmod, misom, p, &pu1);
    1524         217 :   GEN u2 = aux(a3, a4, Tmod, misom, p, &pu2);
    1525         217 :   GEN u3 = aux(a5, a6, Tmod, misom, p, &pu3);
    1526         217 :   pu4 = FpX_mul(pu1,pu2,p);
    1527         217 :   u4 = FpX_chinese_coprime(u1,u2,pu1,pu2,pu4,p);
    1528         217 :   u5 = FpX_chinese_coprime(u4,u3,pu4,pu3,Tp,p);
    1529         217 :   return gerepileupto(av, u5);
    1530             : }
    1531             : static GEN
    1532        5819 : lincomb(GEN A, GEN B, GEN pauto, long j)
    1533             : {
    1534        5819 :   long k = (-j) & 3;
    1535        5819 :   if (j == k) return ZX_mul(ZX_add(A,B), gel(pauto, j+1));
    1536        2951 :   return ZX_add(ZX_mul(A, gel(pauto, j+1)), ZX_mul(B, gel(pauto, k+1)));
    1537             : }
    1538             : /* FIXME: could use the intheadlong technique */
    1539             : static GEN
    1540          35 : s4galoisgen(struct galois_lift *gl)
    1541             : {
    1542          35 :   const long n = 24;
    1543             :   struct galois_testlift gt;
    1544          35 :   pari_sp av, ltop2, ltop = avma;
    1545             :   long i, j;
    1546          35 :   GEN sigma, tau, phi, res, r1,r2,r3,r4, pj, p = gl->p, Q = gl->Q, TQ = gl->TQ;
    1547             :   GEN sg, Tp, Tmod, isom, isominv, misom, Bcoeff, pauto, liftpow, aut;
    1548             : 
    1549          35 :   res = cgetg(3, t_VEC);
    1550          35 :   r1 = cgetg(n+1, t_VECSMALL);
    1551          35 :   r2 = cgetg(n+1, t_VECSMALL);
    1552          35 :   r3 = cgetg(n+1, t_VECSMALL);
    1553          35 :   r4 = cgetg(n+1, t_VECSMALL);
    1554          35 :   gel(res,1)= mkvec4(r1,r2,r3,r4);
    1555          35 :   gel(res,2) = mkvecsmall4(2,2,3,2);
    1556          35 :   ltop2 = avma;
    1557          35 :   sg = identity_perm(6);
    1558          35 :   pj = zero_zv(6);
    1559          35 :   sigma = cgetg(n+1, t_VECSMALL);
    1560          35 :   tau = cgetg(n+1, t_VECSMALL);
    1561          35 :   phi = cgetg(n+1, t_VECSMALL);
    1562          35 :   Tp = FpX_red(gl->T,p);
    1563          35 :   Tmod = gel(FpX_factor(Tp,p), 1);
    1564          35 :   isom    = cgetg(lg(Tmod), t_VEC);
    1565          35 :   isominv = cgetg(lg(Tmod), t_VEC);
    1566          35 :   misom   = cgetg(lg(Tmod), t_MAT);
    1567          35 :   aut = galoisdolift(gl);
    1568          35 :   inittestlift(aut, Tmod, gl, &gt);
    1569          35 :   Bcoeff = gt.bezoutcoeff;
    1570          35 :   pauto = gt.pauto;
    1571         245 :   for (i = 1; i < lg(isom); i++)
    1572             :   {
    1573         210 :     gel(misom,i) = cgetg(lg(Tmod), t_COL);
    1574         210 :     gel(isom,i) = FpX_ffisom(gel(Tmod,1), gel(Tmod,i), p);
    1575         210 :     if (DEBUGLEVEL >= 6)
    1576           0 :       err_printf("S4GaloisConj: Computing isomorphisms %d:%Ps\n", i,
    1577           0 :                  gel(isom,i));
    1578         210 :     gel(isominv,i) = FpXQ_ffisom_inv(gel(isom,i), gel(Tmod,i),p);
    1579             :   }
    1580         245 :   for (i = 1; i < lg(isom); i++)
    1581        1470 :     for (j = 1; j < lg(isom); j++)
    1582        1260 :       gmael(misom,i,j) = FpX_FpXQ_eval(gel(isominv,i),gel(isom,j),
    1583        1260 :                                          gel(Tmod,j),p);
    1584          35 :   liftpow = cgetg(24, t_VEC);
    1585          35 :   av = avma;
    1586          56 :   for (i = 0; i < 3; i++, set_avma(av))
    1587             :   {
    1588             :     pari_sp av1, av2, av3;
    1589             :     GEN u, u1, u2, u3;
    1590             :     long j1, j2, j3;
    1591          56 :     if (i)
    1592             :     {
    1593          21 :       if (i == 1) { lswap(sg[2],sg[3]); }
    1594           0 :       else        { lswap(sg[1],sg[3]); }
    1595             :     }
    1596          56 :     u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[2],sg[3],sg[4],sg[5],sg[6]);
    1597          56 :     s4makelift(u, gl, liftpow);
    1598          56 :     av1 = avma;
    1599         182 :     for (j1 = 0; j1 < 4; j1++, set_avma(av1))
    1600             :     {
    1601         161 :       u1 = lincomb(gel(Bcoeff,sg[5]),gel(Bcoeff,sg[6]), pauto,j1);
    1602         161 :       u1 = FpX_rem(u1, TQ, Q); av2 = avma;
    1603         721 :       for (j2 = 0; j2 < 4; j2++, set_avma(av2))
    1604             :       {
    1605         595 :         u2 = lincomb(gel(Bcoeff,sg[3]),gel(Bcoeff,sg[4]), pauto,j2);
    1606         595 :         u2 = FpX_rem(FpX_add(u1, u2, Q), TQ,Q); av3 = avma;
    1607        2901 :         for (j3 = 0; j3 < 4; j3++, set_avma(av3))
    1608             :         {
    1609        2341 :           u3 = lincomb(gel(Bcoeff,sg[1]),gel(Bcoeff,sg[2]), pauto,j3);
    1610        2341 :           u3 = FpX_rem(FpX_add(u2, u3, Q), TQ,Q);
    1611        2341 :           if (DEBUGLEVEL >= 4)
    1612           0 :             err_printf("S4GaloisConj: Testing %d/3:%d/4:%d/4:%d/4:%Ps\n",
    1613             :                        i,j1,j2,j3, sg);
    1614        2341 :           if (s4test(u3, liftpow, gl, sigma))
    1615             :           {
    1616          35 :             pj[1] = j3;
    1617          35 :             pj[2] = j2;
    1618          35 :             pj[3] = j1; goto suites4;
    1619             :           }
    1620             :         }
    1621             :       }
    1622             :     }
    1623             :   }
    1624           0 :   return gc_NULL(ltop);
    1625          35 : suites4:
    1626          35 :   if (DEBUGLEVEL >= 4) err_printf("S4GaloisConj: sigma=%Ps\n", sigma);
    1627          35 :   if (DEBUGLEVEL >= 4) err_printf("S4GaloisConj: pj=%Ps\n", pj);
    1628          35 :   set_avma(av);
    1629          70 :   for (j = 1; j <= 3; j++)
    1630             :   {
    1631             :     pari_sp av2, av3;
    1632             :     GEN u;
    1633             :     long w, l, z;
    1634          70 :     z = sg[1]; sg[1] = sg[3]; sg[3] = sg[5]; sg[5] = z;
    1635          70 :     z = sg[2]; sg[2] = sg[4]; sg[4] = sg[6]; sg[6] = z;
    1636          70 :     z = pj[1]; pj[1] = pj[2]; pj[2] = pj[3]; pj[3] = z;
    1637         161 :     for (l = 0; l < 2; l++, set_avma(av))
    1638             :     {
    1639         126 :       u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[3],sg[2],sg[4],sg[5],sg[6]);
    1640         126 :       s4makelift(u, gl, liftpow);
    1641         126 :       av2 = avma;
    1642         343 :       for (w = 0; w < 4; w += 2, set_avma(av2))
    1643             :       {
    1644             :         GEN uu;
    1645         252 :         pj[6] = (w + pj[3]) & 3;
    1646         252 :         uu = lincomb(gel(Bcoeff,sg[5]),gel(Bcoeff,sg[6]), pauto, pj[6]);
    1647         252 :         uu = FpX_rem(FpX_red(uu,Q), TQ, Q);
    1648         252 :         av3 = avma;
    1649        1167 :         for (i = 0; i < 4; i++, set_avma(av3))
    1650             :         {
    1651             :           GEN u;
    1652         950 :           pj[4] = i;
    1653         950 :           pj[5] = (i + pj[2] - pj[1]) & 3;
    1654         950 :           if (DEBUGLEVEL >= 4)
    1655           0 :             err_printf("S4GaloisConj: Testing %d/3:%d/2:%d/2:%d/4:%Ps:%Ps\n",
    1656             :                        j-1, w >> 1, l, i, sg, pj);
    1657         950 :           u = ZX_add(lincomb(gel(Bcoeff,sg[1]),gel(Bcoeff,sg[3]), pauto,pj[4]),
    1658         950 :                      lincomb(gel(Bcoeff,sg[2]),gel(Bcoeff,sg[4]), pauto,pj[5]));
    1659         950 :           u = FpX_rem(FpX_add(uu,u,Q), TQ, Q);
    1660         950 :           if (s4test(u, liftpow, gl, tau)) goto suites4_2;
    1661             :         }
    1662             :       }
    1663          91 :       lswap(sg[3], sg[4]);
    1664          91 :       pj[2] = (-pj[2]) & 3;
    1665             :     }
    1666             :   }
    1667           0 :   return gc_NULL(ltop);
    1668          35 : suites4_2:
    1669          35 :   set_avma(av);
    1670             :   {
    1671          35 :     long abc = (pj[1] + pj[2] + pj[3]) & 3;
    1672          35 :     long abcdef = ((abc + pj[4] + pj[5] - pj[6]) & 3) >> 1;
    1673             :     GEN u;
    1674             :     pari_sp av2;
    1675          35 :     u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[4],sg[2],sg[5],sg[3],sg[6]);
    1676          35 :     s4makelift(u, gl, liftpow);
    1677          35 :     av2 = avma;
    1678         190 :     for (j = 0; j < 8; j++)
    1679             :     {
    1680             :       long h, g, i;
    1681         190 :       h = j & 3;
    1682         190 :       g = (abcdef + ((j & 4) >> 1)) & 3;
    1683         190 :       i = (h + abc - g) & 3;
    1684         190 :       u = ZX_add(   lincomb(gel(Bcoeff,sg[1]), gel(Bcoeff,sg[4]), pauto, g),
    1685         190 :                     lincomb(gel(Bcoeff,sg[2]), gel(Bcoeff,sg[5]), pauto, h));
    1686         190 :       u = FpX_add(u, lincomb(gel(Bcoeff,sg[3]), gel(Bcoeff,sg[6]), pauto, i),Q);
    1687         190 :       u = FpX_rem(u, TQ, Q);
    1688         190 :       if (DEBUGLEVEL >= 4)
    1689           0 :         err_printf("S4GaloisConj: Testing %d/8 %d:%d:%d\n", j,g,h,i);
    1690         190 :       if (s4test(u, liftpow, gl, phi)) break;
    1691         155 :       set_avma(av2);
    1692             :     }
    1693             :   }
    1694          35 :   if (j == 8) return gc_NULL(ltop);
    1695         875 :   for (i = 1; i <= n; i++)
    1696             :   {
    1697         840 :     r1[i] = sigma[tau[i]];
    1698         840 :     r2[i] = phi[sigma[tau[phi[i]]]];
    1699         840 :     r3[i] = phi[sigma[i]];
    1700         840 :     r4[i] = sigma[i];
    1701             :   }
    1702          35 :   set_avma(ltop2); return res;
    1703             : }
    1704             : 
    1705             : static GEN
    1706          49 : galoisfindgroups(GEN lo, GEN sg, long f)
    1707             : {
    1708          49 :   pari_sp ltop = avma;
    1709             :   long i, j, k;
    1710          49 :   GEN V = cgetg(lg(lo), t_VEC);
    1711         140 :   for(j=1,i=1; i<lg(lo); i++)
    1712             :   {
    1713          91 :     pari_sp av = avma;
    1714          91 :     GEN loi = gel(lo,i), W = cgetg(lg(loi),t_VECSMALL);
    1715         231 :     for (k=1; k<lg(loi); k++) W[k] = loi[k] % f;
    1716          91 :     W = vecsmall_uniq(W);
    1717          91 :     if (zv_equal(W, sg)) gel(V,j++) = loi;
    1718          91 :     set_avma(av);
    1719             :   }
    1720          49 :   setlg(V,j); return gerepilecopy(ltop,V);
    1721             : }
    1722             : 
    1723             : static GEN
    1724         476 : galoismakepsi(long g, GEN sg, GEN pf)
    1725             : {
    1726         476 :   GEN psi=cgetg(g+1,t_VECSMALL);
    1727             :   long i;
    1728        1176 :   for (i = 1; i < g; i++) psi[i] = sg[pf[i]];
    1729         476 :   psi[g] = sg[1]; return psi;
    1730             : }
    1731             : 
    1732             : static GEN
    1733        2429 : galoisfrobeniuslift_nilp(GEN T, GEN den, GEN L,  GEN Lden,
    1734             :     struct galois_frobenius *gf,  struct galois_borne *gb)
    1735             : {
    1736        2429 :   pari_sp ltop=avma, av2;
    1737             :   struct galois_lift gl;
    1738        2429 :   long i, k, deg = 1, g = lg(gf->Tmod)-1;
    1739        2429 :   GEN F,Fp,Fe, aut, frob, res = cgetg(lg(L), t_VECSMALL);
    1740        2429 :   gf->psi = const_vecsmall(g,1);
    1741        2429 :   av2 = avma;
    1742        2429 :   initlift(T, den, gf->p, L, Lden, gb, &gl);
    1743        2429 :   if (DEBUGLEVEL >= 4)
    1744           0 :     err_printf("GaloisConj: p=%ld e=%ld deg=%ld fp=%ld\n",
    1745             :                             gf->p, gl.e, deg, gf->fp);
    1746        2429 :   aut = galoisdolift(&gl);
    1747        2429 :   if (galoisfrobeniustest(aut,&gl,res))
    1748             :   {
    1749        1890 :     set_avma(av2); gf->deg = gf->fp; return res;
    1750             :   }
    1751             : 
    1752         539 :   F =factoru(gf->fp);
    1753         539 :   Fp = gel(F,1);
    1754         539 :   Fe = gel(F,2);
    1755         539 :   frob = cgetg(lg(L), t_VECSMALL);
    1756        1078 :   for(k = lg(Fp)-1; k>=1; k--)
    1757             :   {
    1758         539 :     pari_sp btop=avma;
    1759         539 :     GEN fres=NULL;
    1760         539 :     long el = gf->fp, dg = 1, dgf = 1, e, pr;
    1761        1078 :     for(e=1; e<=Fe[k]; e++)
    1762             :     {
    1763        1078 :       dg *= Fp[k]; el /= Fp[k];
    1764        1078 :       if (DEBUGLEVEL>=4) err_printf("Trying degre %d.\n",dg);
    1765        1078 :       if (el==1) break;
    1766         637 :       aut = galoisdoliftn(&gl, el);
    1767         637 :       if (!galoisfrobeniustest(aut,&gl,frob))
    1768          98 :         break;
    1769         539 :       dgf = dg; fres = gcopy(frob);
    1770             :     }
    1771         539 :     if (dgf == 1) { set_avma(btop); continue; }
    1772         476 :     pr = deg*dgf;
    1773         476 :     if (deg == 1)
    1774             :     {
    1775       11200 :       for(i=1;i<lg(res);i++) res[i]=fres[i];
    1776             :     }
    1777             :     else
    1778             :     {
    1779           0 :       GEN cp = perm_mul(res,fres);
    1780           0 :       for(i=1;i<lg(res);i++) res[i] = cp[i];
    1781             :     }
    1782         476 :     deg = pr; set_avma(btop);
    1783             :   }
    1784         539 :   if (DEBUGLEVEL>=4 && res) err_printf("Best lift: %d\n",deg);
    1785         539 :   if (deg==1) return gc_NULL(ltop);
    1786             :   else
    1787             :   {
    1788         476 :     set_avma(av2); gf->deg = deg; return res;
    1789             :   }
    1790             : }
    1791             : 
    1792             : static GEN
    1793         658 : galoisfrobeniuslift(GEN T, GEN den, GEN L,  GEN Lden,
    1794             :     struct galois_frobenius *gf,  struct galois_borne *gb)
    1795             : {
    1796         658 :   pari_sp ltop=avma, av2;
    1797             :   struct galois_testlift gt;
    1798             :   struct galois_lift gl;
    1799         658 :   long i, j, k, n = lg(L)-1, deg = 1, g = lg(gf->Tmod)-1;
    1800         658 :   GEN F,Fp,Fe, aut, frob, res = cgetg(lg(L), t_VECSMALL);
    1801         658 :   gf->psi = const_vecsmall(g,1);
    1802         658 :   av2 = avma;
    1803         658 :   initlift(T, den, gf->p, L, Lden, gb, &gl);
    1804         658 :   if (DEBUGLEVEL >= 4)
    1805           0 :     err_printf("GaloisConj: p=%ld e=%ld deg=%ld fp=%ld\n",
    1806             :                             gf->p, gl.e, deg, gf->fp);
    1807         658 :   aut = galoisdolift(&gl);
    1808         658 :   if (galoisfrobeniustest(aut,&gl,res))
    1809             :   {
    1810         182 :     set_avma(av2); gf->deg = gf->fp; return res;
    1811             :   }
    1812         476 :   inittestlift(aut,gf->Tmod, &gl, &gt);
    1813         476 :   gt.C = cgetg(gf->fp+1,t_VEC);
    1814         476 :   gt.Cd= cgetg(gf->fp+1,t_VEC);
    1815        3164 :   for (i = 1; i <= gf->fp; i++) {
    1816        2688 :     gel(gt.C,i)  = zero_zv(gt.g);
    1817        2688 :     gel(gt.Cd,i) = zero_zv(gt.g);
    1818             :   }
    1819             : 
    1820         476 :   F =factoru(gf->fp);
    1821         476 :   Fp = gel(F,1);
    1822         476 :   Fe = gel(F,2);
    1823         476 :   frob = cgetg(lg(L), t_VECSMALL);
    1824        1099 :   for(k=lg(Fp)-1;k>=1;k--)
    1825             :   {
    1826         623 :     pari_sp btop=avma;
    1827         623 :     GEN psi=NULL, fres=NULL, sg = identity_perm(1);
    1828         623 :     long el=gf->fp, dg=1, dgf=1, e, pr;
    1829        1246 :     for(e=1; e<=Fe[k]; e++)
    1830             :     {
    1831             :       GEN lo, pf;
    1832             :       long l;
    1833         672 :       dg *= Fp[k]; el /= Fp[k];
    1834         672 :       if (DEBUGLEVEL>=4) err_printf("Trying degre %d.\n",dg);
    1835         672 :       if (galoisfrobeniustest(gel(gt.pauto,el+1),&gl,frob))
    1836             :       {
    1837         147 :         psi = const_vecsmall(g,1);
    1838         147 :         dgf = dg; fres = gcopy(frob); continue;
    1839             :       }
    1840         525 :       lo = listznstarelts(dg, n / gf->fp);
    1841         525 :       if (e!=1) lo = galoisfindgroups(lo, sg, dgf);
    1842         525 :       if (DEBUGLEVEL>=4) err_printf("Galoisconj:Subgroups list:%Ps\n", lo);
    1843        1057 :       for (l = 1; l < lg(lo); l++)
    1844        1008 :         if (lg(gel(lo,l))>2 && frobeniusliftall(gel(lo,l), el, &pf,&gl,&gt, frob))
    1845             :         {
    1846         476 :           sg  = gcopy(gel(lo,l));
    1847         476 :           psi = galoismakepsi(g,sg,pf);
    1848         476 :           dgf = dg; fres = gcopy(frob); break;
    1849             :         }
    1850         525 :       if (l == lg(lo)) break;
    1851             :     }
    1852         623 :     if (dgf == 1) { set_avma(btop); continue; }
    1853         588 :     pr = deg*dgf;
    1854         588 :     if (deg == 1)
    1855             :     {
    1856        6776 :       for(i=1;i<lg(res);i++) res[i]=fres[i];
    1857        1631 :       for(i=1;i<lg(psi);i++) gf->psi[i]=psi[i];
    1858             :     }
    1859             :     else
    1860             :     {
    1861         112 :       GEN cp = perm_mul(res,fres);
    1862        2422 :       for(i=1;i<lg(res);i++) res[i] = cp[i];
    1863         378 :       for(i=1;i<lg(psi);i++) gf->psi[i] = (dgf*gf->psi[i] + deg*psi[i]) % pr;
    1864             :     }
    1865         588 :     deg = pr; set_avma(btop);
    1866             :   }
    1867        3164 :   for (i = 1; i <= gf->fp; i++)
    1868        8988 :     for (j = 1; j <= gt.g; j++) guncloneNULL(gmael(gt.C,i,j));
    1869         476 :   if (DEBUGLEVEL>=4 && res) err_printf("Best lift: %d\n",deg);
    1870         476 :   if (deg==1) return gc_NULL(ltop);
    1871             :   else
    1872             :   {
    1873             :     /* Normalize result so that psi[g]=1 */
    1874         476 :     long im = Fl_inv(gf->psi[g], deg);
    1875         476 :     GEN cp = perm_pow(res, im);
    1876        6776 :     for(i=1;i<lg(res);i++) res[i] = cp[i];
    1877        1631 :     for(i=1;i<lg(gf->psi);i++) gf->psi[i] = Fl_mul(im, gf->psi[i], deg);
    1878         476 :     set_avma(av2); gf->deg = deg; return res;
    1879             :   }
    1880             : }
    1881             : 
    1882             : /* return NULL if not Galois */
    1883             : static GEN
    1884        3024 : galoisfindfrobenius(GEN T, GEN L, GEN den, GEN bad, struct galois_frobenius *gf,
    1885             :     struct galois_borne *gb, const struct galois_analysis *ga)
    1886             : {
    1887        3024 :   pari_sp ltop = avma, av;
    1888        3024 :   long Try = 0, n = degpol(T), deg, gmask = (ga->group&ga_ext_2)? 3: 1;
    1889        3024 :   GEN frob, Lden = makeLden(L,den,gb);
    1890        3024 :   long is_nilpotent = ga->group&ga_all_nilpotent;
    1891             :   forprime_t S;
    1892             : 
    1893        3024 :   u_forprime_init(&S, ga->p, ULONG_MAX);
    1894        3024 :   av = avma;
    1895        3024 :   deg = gf->deg = ga->deg;
    1896        3087 :   while ((gf->p = u_forprime_next(&S)))
    1897             :   {
    1898             :     pari_sp lbot;
    1899             :     GEN Ti, Tp;
    1900             :     long nb, d;
    1901        3087 :     set_avma(av);
    1902        3087 :     Tp = ZX_to_Flx(T, gf->p);
    1903        3087 :     if (!Flx_is_squarefree(Tp, gf->p)) continue;
    1904        3087 :     if (bad && dvdiu(bad, gf->p)) continue;
    1905        3087 :     Ti = gel(Flx_factor(Tp, gf->p), 1);
    1906        3087 :     nb = lg(Ti)-1; d = degpol(gel(Ti,1));
    1907        3087 :     if (nb > 1 && degpol(gel(Ti,nb)) != d) return gc_NULL(ltop);
    1908        3087 :     if (((gmask&1)==0 || d % deg) && ((gmask&2)==0 || odd(d))) continue;
    1909        3087 :     if (DEBUGLEVEL >= 1) err_printf("GaloisConj: Trying p=%ld\n", gf->p);
    1910        3087 :     FlxV_to_ZXV_inplace(Ti);
    1911        3087 :     gf->fp = d;
    1912        3087 :     gf->Tmod = Ti; lbot = avma;
    1913        3087 :     if (is_nilpotent)
    1914        2429 :       frob = galoisfrobeniuslift_nilp(T, den, L, Lden, gf, gb);
    1915             :     else
    1916         658 :       frob = galoisfrobeniuslift(T, den, L, Lden, gf, gb);
    1917        3087 :     if (frob)
    1918             :     {
    1919             :       GEN *gptr[3];
    1920        3024 :       gf->Tmod = gcopy(Ti);
    1921        3024 :       gptr[0]=&gf->Tmod; gptr[1]=&gf->psi; gptr[2]=&frob;
    1922        3024 :       gerepilemanysp(ltop,lbot,gptr,3); return frob;
    1923             :     }
    1924          63 :     if (is_nilpotent) continue;
    1925           0 :     if ((ga->group&ga_all_normal) && d % deg == 0) gmask &= ~1;
    1926             :     /* The first prime degree is always divisible by deg, so we don't
    1927             :      * have to worry about ext_2 being used before regular supersolvable*/
    1928           0 :     if (!gmask) return gc_NULL(ltop);
    1929           0 :     if ((ga->group&ga_non_wss) && ++Try > ((3*n)>>1))
    1930             :     {
    1931           0 :       pari_warn(warner,"Galois group probably not weakly super solvable");
    1932           0 :       return NULL;
    1933             :     }
    1934             :   }
    1935           0 :   if (!gf->p) pari_err_OVERFLOW("galoisfindfrobenius [ran out of primes]");
    1936           0 :   return NULL;
    1937             : }
    1938             : 
    1939             : /* compute g such that tau(Pmod[#])= tau(Pmod[g]) */
    1940             : 
    1941             : static long
    1942        3171 : get_image(GEN tau, GEN P, GEN Pmod, GEN p)
    1943             : {
    1944        3171 :   pari_sp av = avma;
    1945        3171 :   long g, gp = lg(Pmod)-1;
    1946        3171 :   tau = RgX_to_FpX(tau, p);
    1947        3171 :   tau = FpX_FpXQ_eval(gel(Pmod, gp), tau, P, p);
    1948        3171 :   tau = FpX_normalize(FpX_gcd(P, tau, p), p);
    1949        6601 :   for (g = 1; g <= gp; g++)
    1950        6601 :     if (ZX_equal(tau, gel(Pmod,g))) return gc_long(av,g);
    1951           0 :   return gc_long(av,0);
    1952             : }
    1953             : 
    1954             : static GEN
    1955        4487 : gg_get_std(GEN G)
    1956             : {
    1957        4487 :   return !G ? NULL: lg(G)==3 ? G: mkvec2(gel(G,1),gmael(G,5,1));
    1958             : }
    1959             : 
    1960             : static GEN galoisgen(GEN T, GEN L, GEN M, GEN den, GEN bad, struct galois_borne *gb,
    1961             :           const struct galois_analysis *ga);
    1962             : 
    1963             : static GEN
    1964        2163 : galoisgenfixedfield(GEN Tp, GEN Pmod, GEN PL, GEN P, GEN ip, GEN bad, struct galois_borne *gb)
    1965             : {
    1966             :   GEN  Pden, PM;
    1967             :   GEN  tau, PG, Pg;
    1968             :   long g, lP;
    1969        2163 :   long x = varn(Tp);
    1970        2163 :   GEN Pp = FpX_red(P, ip);
    1971        2163 :   if (DEBUGLEVEL>=6)
    1972           0 :     err_printf("GaloisConj: Fixed field %Ps\n",P);
    1973        2163 :   if (degpol(P)==2 && !bad)
    1974             :   {
    1975        1393 :     PG=cgetg(3,t_VEC);
    1976        1393 :     gel(PG,1) = mkvec( mkvecsmall2(2,1) );
    1977        1393 :     gel(PG,2) = mkvecsmall(2);
    1978        1393 :     tau = deg1pol_shallow(gen_m1, negi(gel(P,3)), x);
    1979        1393 :     g = get_image(tau, Pp, Pmod, ip);
    1980        1393 :     if (!g) return NULL;
    1981        1393 :     Pg = mkvecsmall(g);
    1982             :   }
    1983             :   else
    1984             :   {
    1985             :     struct galois_analysis Pga;
    1986             :     struct galois_borne Pgb;
    1987             :     GEN mod, mod2;
    1988             :     long j;
    1989         777 :     if (!galoisanalysis(P, &Pga, 0, NULL)) return NULL;
    1990         756 :     if (bad) Pga.group &= ~ga_easy;
    1991         756 :     Pgb.l = gb->l;
    1992         756 :     Pden = galoisborne(P, NULL, &Pgb, degpol(P));
    1993             : 
    1994         756 :     if (Pgb.valabs > gb->valabs)
    1995             :     {
    1996         118 :       if (DEBUGLEVEL>=4)
    1997           0 :         err_printf("GaloisConj: increase prec of p-adic roots of %ld.\n"
    1998           0 :             ,Pgb.valabs-gb->valabs);
    1999         118 :       PL = ZpX_liftroots(P,PL,gb->l,Pgb.valabs);
    2000             :     }
    2001         638 :     else if (Pgb.valabs < gb->valabs)
    2002         574 :       PL = FpC_red(PL, Pgb.ladicabs);
    2003         756 :     PM = FpV_invVandermonde(PL, Pden, Pgb.ladicabs);
    2004         756 :     PG = galoisgen(P, PL, PM, Pden, bad ? lcmii(Pgb.dis, bad): NULL, &Pgb, &Pga);
    2005         756 :     if (!PG) return NULL;
    2006         749 :     lP = lg(gel(PG,1));
    2007         749 :     mod = Pgb.ladicabs; mod2 = shifti(mod, -1);
    2008         749 :     Pg = cgetg(lP, t_VECSMALL);
    2009        2527 :     for (j = 1; j < lP; j++)
    2010             :     {
    2011        1778 :       pari_sp btop=avma;
    2012        1778 :       tau = permtopol(gmael(PG,1,j), PL, PM, Pden, mod, mod2, x);
    2013        1778 :       g = get_image(tau, Pp, Pmod, ip);
    2014        1778 :       if (!g) return NULL;
    2015        1778 :       Pg[j] = g;
    2016        1778 :       set_avma(btop);
    2017             :     }
    2018             :   }
    2019        2142 :   return mkvec2(PG,Pg);
    2020             : }
    2021             : 
    2022             : static GEN
    2023        2163 : galoisgenfixedfield0(GEN O, GEN L, GEN sigma, GEN T, GEN bad, GEN *pt_V,
    2024             :                      struct galois_frobenius *gf, struct galois_borne *gb)
    2025             : {
    2026        2163 :   pari_sp btop = avma;
    2027        2163 :   long vT = varn(T);
    2028        2163 :   GEN mod = gb->ladicabs, mod2 = shifti(gb->ladicabs,-1);
    2029             :   GEN OL, sym, P, PL, p, Tp, Sp, Pmod, PG;
    2030        2163 :   OL = fixedfieldorbits(O,L);
    2031        2163 :   sym  = fixedfieldsympol(OL, itou(gb->l));
    2032        2163 :   PL = sympol_eval(sym, OL, mod);
    2033        2163 :   P = FpX_center_i(FpV_roots_to_pol(PL, mod, vT), mod, mod2);
    2034        2163 :   if (!FpX_is_squarefree(P,utoipos(gf->p)))
    2035             :   {
    2036          70 :     GEN badp = lcmii(bad? bad: gb->dis, ZX_disc(P));
    2037          70 :     gf->p  = findpsi(badp, gf->p, T, sigma, gf->deg, &gf->Tmod, &gf->psi);
    2038             :   }
    2039        2163 :   p  = utoipos(gf->p);
    2040        2163 :   Tp = FpX_red(T,p);
    2041        2163 :   Sp = sympol_aut_evalmod(sym, gf->deg, sigma, Tp, p);
    2042        2163 :   Pmod = fixedfieldfactmod(Sp, p, gf->Tmod);
    2043        2163 :   PG = galoisgenfixedfield(Tp, Pmod, PL, P, p, bad, gb);
    2044        2163 :   if (PG == NULL) return NULL;
    2045        2142 :   if (DEBUGLEVEL >= 4)
    2046           0 :     err_printf("GaloisConj: Back to Earth:%Ps\n", gg_get_std(gel(PG,1)));
    2047        2142 :   if (pt_V) *pt_V = mkvec3(sym, PL, P);
    2048        2142 :   gerepileall(btop, pt_V ? 4: 3, &gf->Tmod, &gf->psi, &PG, pt_V);
    2049        2142 :   return PG;
    2050             : }
    2051             : 
    2052             : /* Let sigma^m=1, tau*sigma*tau^-1=sigma^s. Return n = sum_{0<=k<e,0} s^k mod m
    2053             :  * so that (sigma*tau)^e = sigma^n*tau^e. N.B. n*(1-s) = 1-s^e mod m,
    2054             :  * unfortunately (1-s) may not invertible mod m */
    2055             : static long
    2056        7056 : stpow(long s, long e, long m)
    2057             : {
    2058        7056 :   long i, n = 1;
    2059       11088 :   for (i = 1; i < e; i++) n = (1 + n * s) % m;
    2060        7056 :   return n;
    2061             : }
    2062             : 
    2063             : static GEN
    2064        3171 : wpow(long s, long m, long e, long n)
    2065             : {
    2066        3171 :   GEN   w = cgetg(n+1,t_VECSMALL);
    2067        3171 :   long si = s;
    2068             :   long i;
    2069        3171 :   w[1] = 1;
    2070        3528 :   for(i=2; i<=n; i++) w[i] = w[i-1]*e;
    2071        6699 :   for(i=n; i>=1; i--)
    2072             :   {
    2073        3528 :     si = Fl_powu(si,e,m);
    2074        3528 :     w[i] = Fl_mul(s-1, stpow(si, w[i], m), m);
    2075             :   }
    2076        3171 :   return w;
    2077             : }
    2078             : 
    2079             : static GEN
    2080        3171 : galoisgenliftauto(GEN O, GEN gj, long s, long n, struct galois_test *td)
    2081             : {
    2082        3171 :   pari_sp av = avma;
    2083             :   long sr, k;
    2084        3171 :   long deg = lg(gel(O,1))-1;
    2085        3171 :   GEN  X  = cgetg(lg(O), t_VECSMALL);
    2086        3171 :   GEN  oX = cgetg(lg(O), t_VECSMALL);
    2087        3171 :   GEN  B  = perm_cycles(gj);
    2088        3171 :   long oj = lg(gel(B,1)) - 1;
    2089        3171 :   GEN  F  = factoru(oj);
    2090        3171 :   GEN  Fp = gel(F,1);
    2091        3171 :   GEN  Fe = gel(F,2);
    2092        3171 :   GEN  pf = identity_perm(n);
    2093        3171 :   if (DEBUGLEVEL >= 6)
    2094           0 :     err_printf("GaloisConj: %Ps of relative order %d\n", gj, oj);
    2095        6321 :   for (k=lg(Fp)-1; k>=1; k--)
    2096             :   {
    2097        3171 :     long f, dg = 1, el = oj, osel = 1, a = 0;
    2098        3171 :     long p  = Fp[k], e  = Fe[k], op = oj / upowuu(p,e);
    2099             :     long i;
    2100        3171 :     GEN  pf1 = NULL, w, wg, Be = cgetg(e+1,t_VEC);
    2101        3171 :     gel(Be,e) = cyc_pow(B, op);
    2102        3528 :     for(i=e-1; i>=1; i--) gel(Be,i) = cyc_pow(gel(Be,i+1), p);
    2103        3171 :     w = wpow(Fl_powu(s,op,deg),deg,p,e);
    2104        3171 :     wg = cgetg(e+2,t_VECSMALL);
    2105        3171 :     wg[e+1] = deg;
    2106        6699 :     for (i=e; i>=1; i--) wg[i] = ugcd(wg[i+1], w[i]);
    2107       24045 :     for (i=1; i<lg(O); i++) oX[i] = 0;
    2108        6678 :     for (f=1; f<=e; f++)
    2109             :     {
    2110             :       long sel, t;
    2111        3528 :       GEN Bel = gel(Be,f);
    2112        3528 :       dg *= p; el /= p;
    2113        3528 :       sel = Fl_powu(s,el,deg);
    2114        3528 :       if (DEBUGLEVEL >= 6) err_printf("GaloisConj: B=%Ps\n", Bel);
    2115        3528 :       sr  = ugcd(stpow(sel,p,deg),deg);
    2116        3528 :       if (DEBUGLEVEL >= 6)
    2117           0 :         err_printf("GaloisConj: exp %d: s=%ld [%ld] a=%ld w=%ld wg=%ld sr=%ld\n",
    2118           0 :             dg, sel, deg, a, w[f], wg[f+1], sr);
    2119        4424 :       for (t = 0; t < sr; t++)
    2120        4403 :         if ((a+t*w[f])%wg[f+1]==0)
    2121             :         {
    2122             :           long i, j, k, st;
    2123       41055 :           for (i = 1; i < lg(X); i++) X[i] = 0;
    2124       20181 :           for (i = 0; i < lg(X)-1; i+=dg)
    2125       34111 :             for (j = 1, k = p, st = t; k <= dg; j++, k += p)
    2126             :             {
    2127       18263 :               X[k+i] = (oX[j+i] + st)%deg;
    2128       18263 :               st = (t + st*osel)%deg;
    2129             :             }
    2130        4333 :           pf1 = testpermutation(O, Bel, X, sel, p, sr, td);
    2131        4333 :           if (pf1) break;
    2132             :         }
    2133        3528 :       if (!pf1) return NULL;
    2134       30583 :       for (i=1; i<lg(O); i++) oX[i] = X[i];
    2135        3507 :       osel = sel; a = (a+t*w[f])%deg;
    2136             :     }
    2137        3150 :     pf = perm_mul(pf, perm_pow(pf1, el));
    2138             :   }
    2139        3150 :   return gerepileuptoleaf(av, pf);
    2140             : }
    2141             : 
    2142             : static GEN
    2143           0 : FlxV_Flx_gcd(GEN x, GEN T, ulong p)
    2144           0 : { pari_APPLY_same(Flx_normalize(Flx_gcd(gel(x,i),T,p),p)) }
    2145             : 
    2146             : static GEN
    2147           0 : Flx_FlxV_minpolymod(GEN y, GEN x, ulong p)
    2148           0 : { pari_APPLY_same(Flxq_minpoly(Flx_rem(y, gel(x,i), p), gel(x,i), p)) }
    2149             : 
    2150             : static GEN
    2151           0 : FlxV_minpolymod(GEN x, GEN y, ulong p)
    2152           0 : { pari_APPLY_same(Flx_FlxV_minpolymod(gel(x,i), y, p)) }
    2153             : 
    2154             : static GEN
    2155           0 : factperm(GEN x)
    2156             : {
    2157           0 :   pari_APPLY_same(gen_indexsort(gel(x,i), (void*)cmp_Flx, cmp_nodata))
    2158             : }
    2159             : 
    2160             : /* compute (prod p_i^e_i)(1) */
    2161             : 
    2162             : static long
    2163           0 : permprodeval(GEN p, GEN e, long s)
    2164             : {
    2165           0 :   long i, j, l = lg(p);
    2166           0 :   for (i=l-1; i>=1; i--)
    2167             :   {
    2168           0 :     GEN pi = gel(p,i);
    2169           0 :     long ei = uel(e,i);
    2170           0 :     for(j = 1; j <= ei; j++)
    2171           0 :       s = uel(pi, s);
    2172             :   }
    2173           0 :   return s;
    2174             : }
    2175             : 
    2176             : static GEN
    2177           0 : pc_to_perm(GEN pc, GEN gen, long n)
    2178             : {
    2179           0 :   long i, l = lg(pc);
    2180           0 :   GEN s = identity_perm(n);
    2181           0 :   for (i=1; i<l; i++)
    2182           0 :     s = perm_mul(gel(gen,pc[i]),s);
    2183           0 :   return s;
    2184             : }
    2185             : 
    2186             : static GEN
    2187           0 : genorbit(GEN ordH, GEN permfact_Hp, long fr, long n, long k, long j)
    2188             : {
    2189           0 :   pari_sp av = avma;
    2190           0 :   long l = lg(gel(permfact_Hp,1))-1, no = 1, b, i;
    2191           0 :   GEN W = zero_zv(l);
    2192           0 :   GEN orb = cgetg(l+1, t_VECSMALL);
    2193           0 :   GEN gen = cgetg(l+1, t_VEC);
    2194           0 :   GEN E = cgetg(k+1, t_VECSMALL);
    2195           0 :   for(b = 0; b < n; b++)
    2196             :   {
    2197           0 :     long bb = b, s;
    2198           0 :     for(i = 1; i <= k; i++)
    2199             :     {
    2200           0 :       uel(E,i) = bb % uel(ordH,i);
    2201           0 :       bb /= uel(ordH,i);
    2202             :     }
    2203           0 :     if (E[j]) continue;
    2204           0 :     s = permprodeval(permfact_Hp, E, fr);
    2205           0 :     if (s>lg(W)-1) pari_err_BUG("W1");
    2206           0 :     if (W[s]) continue;
    2207           0 :     W[s] = 1;
    2208           0 :     if (no > l) pari_err_BUG("genorbit");
    2209           0 :     uel(orb,no) = s;
    2210           0 :     gel(gen,no) = zv_copy(E);
    2211           0 :     no++;
    2212             :   }
    2213           0 :   if(no<l) pari_err_BUG("genorbit");
    2214           0 :   return gerepilecopy(av, mkvec2(orb,gen));
    2215             : }
    2216             : 
    2217           0 : INLINE GEN br_get(GEN br, long i, long j) { return gmael(br,j,i-j); }
    2218           0 : static GEN pcgrp_get_ord(GEN G) { return gel(G,1); }
    2219           0 : static GEN pcgrp_get_pow(GEN G) { return gel(G,2); }
    2220           0 : static GEN pcgrp_get_br(GEN G)  { return gel(G,3); }
    2221             : 
    2222             : static GEN
    2223         861 : cyclic_pc(long n)
    2224             : {
    2225         861 :   return mkvec3(mkvecsmall(n),mkvec(cgetg(1,t_VECSMALL)), mkvec(cgetg(1,t_VEC)));
    2226             : }
    2227             : 
    2228             : static GEN
    2229           0 : pc_normalize(GEN g, GEN G)
    2230             : {
    2231           0 :   long i, l = lg(g)-1, o = 1;
    2232           0 :   GEN ord = pcgrp_get_ord(G), pw = pcgrp_get_pow(G), br = pcgrp_get_br(G);
    2233           0 :   for (i = 1; i < l; i++)
    2234             :   {
    2235           0 :     if (g[i] == g[i+1])
    2236             :     {
    2237           0 :       if (++o == ord[g[i]])
    2238             :       {
    2239           0 :         GEN v = vecsmall_concat(vecslice(g,1,i-o+1),gel(pw,g[i]));
    2240           0 :         GEN w = vecsmall_concat(v,vecslice(g,i+2,l));
    2241           0 :         return pc_normalize(w, G);
    2242             :       }
    2243             :     }
    2244           0 :     else if (g[i] > g[i+1])
    2245             :     {
    2246           0 :       GEN v = vecsmall_concat(vecslice(g,1,i-1), br_get(br,g[i],g[i+1]));
    2247           0 :       GEN w = vecsmall_concat(mkvecsmall2(g[i+1],g[i]),vecslice(g,i+2,l));
    2248           0 :       v = vecsmall_concat(v, w);
    2249           0 :       return pc_normalize(v, G);
    2250             :     }
    2251           0 :     else o = 1;
    2252             :   }
    2253           0 :   return g;
    2254             : }
    2255             : 
    2256             : static GEN
    2257           0 : pc_inv(GEN g, GEN G)
    2258             : {
    2259           0 :   long i, l = lg(g);
    2260           0 :   GEN ord = pcgrp_get_ord(G), pw  = pcgrp_get_pow(G);
    2261           0 :   GEN v = cgetg(l, t_VEC);
    2262           0 :   if (l==1) return v;
    2263           0 :   for(i = 1; i < l; i++)
    2264             :   {
    2265           0 :     ulong gi = uel(g,i);
    2266           0 :     gel(v,l-i) = vecsmall_concat(pc_inv(gel(pw, gi), G),
    2267           0 :                                  const_vecsmall(uel(ord,gi)-1,gi));
    2268             :   }
    2269           0 :   return pc_normalize(shallowconcat1(v), G);
    2270             : }
    2271             : 
    2272             : static GEN
    2273           0 : pc_mul(GEN g, GEN h, GEN G)
    2274             : {
    2275           0 :   return pc_normalize(vecsmall_concat(g,h), G);
    2276             : }
    2277             : 
    2278             : static GEN
    2279           0 : pc_bracket(GEN g, GEN h, GEN G)
    2280             : {
    2281           0 :   GEN gh = pc_mul(g, h, G);
    2282           0 :   GEN hg = pc_mul(h, g, G);
    2283           0 :   long i, l1 = lg(gh), l2 = lg(hg), lm = minss(l1,l2);
    2284           0 :   for (i = 1; i < lm; i++)
    2285           0 :     if (gh[l1-i] != hg[l2-i]) break;
    2286           0 :   return pc_mul(vecsmall_shorten(gh,l1-i), pc_inv(vecsmall_shorten(hg,l2-i), G), G);
    2287             : }
    2288             : 
    2289             : static GEN
    2290           0 : pc_exp(GEN v)
    2291             : {
    2292           0 :   long i, l = lg(v);
    2293           0 :   GEN w = cgetg(l, t_VEC);
    2294           0 :   if (l==1) return w;
    2295           0 :   for (i = 1; i < l; i++)
    2296           0 :     gel(w,i) = const_vecsmall(v[i], i+1);
    2297           0 :   return shallowconcat1(w);
    2298             : }
    2299             : static GEN
    2300           0 : vecsmall_increase(GEN x)
    2301           0 : { pari_APPLY_ulong(x[i]+1) }
    2302             : 
    2303             : static GEN
    2304           0 : vecvecsmall_increase(GEN x)
    2305           0 : { pari_APPLY_same(vecsmall_increase(gel(x,i))) }
    2306             : 
    2307             : static GEN
    2308           0 : pcgrp_lift(GEN G, long deg)
    2309             : {
    2310           0 :   GEN ord = pcgrp_get_ord(G), pw  = pcgrp_get_pow(G), br = pcgrp_get_br(G);
    2311           0 :   long i, l = lg(br);
    2312           0 :   GEN Ord = vecsmall_prepend(ord, deg);
    2313           0 :   GEN Pw = vec_prepend(vecvecsmall_increase(pw), cgetg(1,t_VECSMALL));
    2314           0 :   GEN Br = cgetg(l+1, t_VEC);
    2315           0 :   gel(Br,1) = const_vec(l-1, cgetg(1, t_VECSMALL));
    2316           0 :   for (i = 1; i < l; i++)
    2317           0 :     gel(Br,i+1) = vecvecsmall_increase(gel(br, i));
    2318           0 :   return mkvec3(Ord, Pw, Br);
    2319             : }
    2320             : 
    2321             : static GEN
    2322           0 : brl_add(GEN x, GEN a)
    2323             : {
    2324           0 :   pari_APPLY_same(vecsmall_concat(const_vecsmall(uel(a,i),1),gel(x,i)))
    2325             : }
    2326             : 
    2327             : static void
    2328           0 : pcgrp_insert(GEN G, long j, GEN a)
    2329             : {
    2330           0 :   GEN pw  = pcgrp_get_pow(G), br = pcgrp_get_br(G);
    2331           0 :   gel(pw,j) = vecsmall_concat(gel(a,1),gel(pw, j));
    2332           0 :   gel(br,j) = brl_add(gel(br, j), gel(a,2));
    2333           0 : }
    2334             : 
    2335             : static long
    2336           0 : getfr(GEN f, GEN h)
    2337             : {
    2338           0 :   long i, l = lg(f);
    2339           0 :   for (i = 1; i < l; i++)
    2340           0 :     if (zv_equal(gel(f,i), h)) return i;
    2341           0 :   pari_err_BUG("galoisinit");
    2342           0 :   return 0;
    2343             : }
    2344             : 
    2345             : static long
    2346           0 : get_pow(GEN pf, long o, GEN pw, GEN gen)
    2347             : {
    2348           0 :   long i, n  = lg(pf)-1;
    2349           0 :   GEN p1 = perm_pow(pf, o);
    2350           0 :   GEN p2 = pc_to_perm(pw, gen, n);
    2351           0 :   for(i = 0; ; i++)
    2352             :   {
    2353           0 :     if (zv_equal(p1, p2)) break;
    2354           0 :     p2 = perm_mul(gel(gen,1), p2);
    2355             :   }
    2356           0 :   return i;
    2357             : }
    2358             : 
    2359             : struct galois_perm
    2360             : {
    2361             :   GEN L;
    2362             :   GEN M;
    2363             :   GEN den;
    2364             :   GEN mod, mod2;
    2365             :   long x;
    2366             :   GEN cache;
    2367             : };
    2368             : 
    2369             : static void
    2370           0 : galoisperm_init(struct galois_perm *gp, GEN L, GEN M, GEN den, GEN mod, GEN mod2, long x)
    2371             : {
    2372           0 :   gp->L = L;
    2373           0 :   gp->M = M;
    2374           0 :   gp->den = den;
    2375           0 :   gp->mod = mod;
    2376           0 :   gp->mod2 = mod2;
    2377           0 :   gp->x = x;
    2378           0 :   gp->cache = zerovec(lg(L)-1);
    2379           0 : }
    2380             : 
    2381             : static void
    2382           0 : galoisperm_free(struct galois_perm *gp)
    2383             : {
    2384           0 :   long i, l = lg(gp->cache);
    2385           0 :   for (i=1; i<l; i++)
    2386           0 :     if (!isintzero(gel(gp->cache,i)))
    2387           0 :       gunclone(gel(gp->cache,i));
    2388           0 : }
    2389             : 
    2390             : static GEN
    2391           0 : permtoaut(GEN p, struct galois_perm *gp)
    2392             : {
    2393           0 :   pari_sp av = avma;
    2394           0 :   if (isintzero(gel(gp->cache,p[1])))
    2395             :   {
    2396           0 :     GEN pol = permtopol(p, gp->L, gp->M, gp->den, gp->mod, gp->mod2, gp->x);
    2397           0 :     gel(gp->cache,p[1]) = gclone(pol);
    2398             :   }
    2399           0 :   set_avma(av);
    2400           0 :   return gel(gp->cache,p[1]);
    2401             : }
    2402             : 
    2403             : static GEN
    2404           0 : pc_evalcache(GEN W, GEN u, GEN ss, GEN genG, GEN T, GEN p, struct galois_perm *gp)
    2405             : {
    2406             :   GEN sp, v;
    2407             :   long ns;
    2408           0 :   if (lg(ss) == 1) return gel(u,2);
    2409           0 :   sp = pc_to_perm(ss, genG, degpol(T));
    2410           0 :   ns = sp[1];
    2411           0 :   if (!isintzero(gel(W,ns))) return gel(W,ns);
    2412           0 :   v = RgX_to_FpX(permtoaut(sp, gp), p);
    2413           0 :   gel(W,ns) = FpX_FpXQV_eval(v, u, T, p);
    2414           0 :   return gel(W,ns);
    2415             : }
    2416             : 
    2417             : static ulong
    2418           0 : findp(GEN D, GEN P, GEN S, long o, GEN *Tmod)
    2419             : {
    2420             :   forprime_t iter;
    2421             :   ulong p;
    2422           0 :   long n = degpol(P);
    2423           0 :   u_forprime_init(&iter, n*maxss(expu(n)-3, 2), ULONG_MAX);
    2424           0 :   while ((p = u_forprime_next(&iter)))
    2425             :   {
    2426             :     GEN F, F1, Sp;
    2427           0 :     if (smodis(D, p) == 0)
    2428           0 :       continue;
    2429           0 :     F = gel(Flx_factor(ZX_to_Flx(P, p), p), 1);
    2430           0 :     F1 = gel(F,1);
    2431           0 :     if (degpol(F1) != o)
    2432           0 :       continue;
    2433           0 :     Sp = RgX_to_Flx(S, p);
    2434           0 :     if (gequal(Flx_rem(Sp, F1, p), Flx_Frobenius(F1, p)))
    2435             :     {
    2436           0 :       *Tmod = FlxV_to_ZXV(F);
    2437           0 :       return p;
    2438             :     }
    2439             :   }
    2440           0 :   return 0;
    2441             : }
    2442             : 
    2443             : static GEN
    2444           0 : nilp_froblift(GEN genG, GEN autH, long j, GEN pcgrp,
    2445             :   GEN idp, GEN incl, GEN H, struct galois_lift *gl, struct galois_perm *gp)
    2446             : {
    2447           0 :   pari_sp av = avma;
    2448           0 :   GEN T = gl->T, p = gl->p, pe = gl->Q;
    2449           0 :   ulong pp = itou(p);
    2450           0 :   long e   = gl->e;
    2451           0 :   GEN pf   = cgetg(lg(gl->L), t_VECSMALL);
    2452           0 :   GEN Tp   = ZX_to_Flx(T, pp);
    2453           0 :   GEN Hp   = ZX_to_Flx(H, pp);
    2454           0 :   GEN ord = pcgrp_get_ord(pcgrp);
    2455           0 :   GEN pcp = gel(pcgrp_get_pow(pcgrp),j+1);
    2456           0 :   long o  = uel(ord,1);
    2457           0 :   GEN ordH = vecslice(ord,2,lg(ord)-1);
    2458           0 :   long n = zv_prod(ordH), k = lg(ordH)-1, l = k-j, m = upowuu(o, l), v = varn(T);
    2459           0 :   GEN factTp = gel(Flx_factor(Tp, pp), 1);
    2460           0 :   long fp = degpol(gel(factTp, 1));
    2461           0 :   GEN frobp = Flxq_autpow(Flx_Frobenius(Tp, pp), fp-1, Tp, pp);
    2462           0 :   GEN frob = ZpX_ZpXQ_liftroot(T, Flx_to_ZX(frobp), T, p, e);
    2463           0 :   if (galoisfrobeniustest(frob, gl, pf))
    2464             :   {
    2465           0 :     GEN pfi = perm_inv(pf);
    2466           0 :     long d = get_pow(pfi, uel(ord,j+1), pcp, genG);
    2467           0 :     return mkvec3(pfi, mkvec2(const_vecsmall(d,1),zero_zv(l+1)), gel(factTp, 1));
    2468             :   }
    2469             :   else
    2470             :   {
    2471           0 :     GEN frobG = FpXQ_powers(frob, usqrt(degpol(T)), T, pe);
    2472           0 :     GEN autHp = RgXV_to_FlxV(autH,pp);
    2473           0 :     GEN inclp = RgX_to_Flx(incl,pp);
    2474           0 :     GEN factHp = gel(Flx_factor(Hp, pp),1);
    2475           0 :     long fr = getfr(factHp, idp);
    2476           0 :     GEN minHp  = FlxV_minpolymod(autHp, factHp, pp);
    2477           0 :     GEN permfact_Hp = factperm(minHp);
    2478           0 :     GEN permfact_Gp = FlxV_Flx_gcd(FlxC_Flxq_eval(factHp, inclp, Tp, pp), Tp, pp);
    2479           0 :     GEN bezout_Gpe = bezout_lift_fact(T, FlxV_to_ZXV(permfact_Gp), p, e);
    2480           0 :     GEN id = gmael(Flx_factor(gel(permfact_Gp, fr),pp),1,1);
    2481           0 :     GEN orbgen = genorbit(ordH, permfact_Hp, fr, n, k, j);
    2482           0 :     GEN orb = gel(orbgen,1), gen = gel(orbgen,2);
    2483           0 :     long nborb = lg(orb)-1;
    2484           0 :     GEN A = cgetg(l+1, t_VECSMALL);
    2485           0 :     GEN W = zerovec(lg(gl->L)-1);
    2486           0 :     GEN br = pcgrp_get_br(pcgrp), brj = gcopy(gel(br, j+1));
    2487           0 :     GEN U = cgetg(nborb+1, t_VEC);
    2488             :     long a, b, i;
    2489           0 :     for(a = 0; a < m; a++)
    2490             :     {
    2491             :       pari_sp av2;
    2492           0 :       GEN B = pol_0(v);
    2493           0 :       long aa = a;
    2494           0 :       for(i = 1; i <= l; i++)
    2495             :       {
    2496           0 :         uel(A,i) = aa % o;
    2497           0 :         aa /= o;
    2498             :       }
    2499           0 :       gel(br,j+1) = brl_add(brj, A);
    2500           0 :       for(b = 1; b <= nborb; b++)
    2501             :       {
    2502           0 :         GEN br = pc_bracket(pc_exp(gel(gen,b)), mkvecsmall(j+1), pcgrp);
    2503           0 :         gel(U, b) = pc_evalcache(W, frobG, br, genG, T, pe, gp);
    2504             :       }
    2505           0 :       av2 = avma;
    2506           0 :       for(b = 1; b <= nborb; b++)
    2507             :       {
    2508           0 :         long s = permprodeval(permfact_Hp, gel(gen,b), fr);
    2509           0 :         B = FpX_add(B, FpXQ_mul(gel(U, b), gel(bezout_Gpe,s), T, pe), pe);
    2510             :       }
    2511           0 :       if (galoisfrobeniustest(B, gl, pf))
    2512             :       {
    2513           0 :         GEN pfi = perm_inv(pf);
    2514           0 :         long d = get_pow(pfi, uel(ord,j+1), pcp, genG);
    2515           0 :         gel(br,j+1) = brj;
    2516           0 :         return gerepilecopy(av,mkvec3(pfi,mkvec2(const_vecsmall(d,1),A),id));
    2517             :       }
    2518           0 :       set_avma(av2);
    2519             :     }
    2520           0 :     return gc_NULL(av);
    2521             :   }
    2522             : }
    2523             : 
    2524             : static GEN
    2525           0 : galoisgenlift_nilp(GEN PG, GEN O, GEN V, GEN T, GEN frob, GEN sigma,
    2526             :   struct galois_borne *gb, struct galois_frobenius *gf, struct galois_perm *gp)
    2527             : {
    2528           0 :   long j, n = degpol(T), deg = gf->deg;
    2529           0 :   ulong p = gf->p;
    2530           0 :   GEN L = gp->L, M =  gp->M, den = gp->den;
    2531           0 :   GEN S = fixedfieldinclusion(O, gel(V,2));
    2532           0 :   GEN incl = vectopol(S, M, den, gb->ladicabs, shifti(gb->ladicabs,-1), varn(T));
    2533           0 :   GEN H = gel(V,3);
    2534           0 :   GEN PG1 = gmael(PG, 1, 1);
    2535           0 :   GEN PG2 = gmael(PG, 1, 2);
    2536           0 :   GEN PG3 = gmael(PG, 1, 3);
    2537           0 :   GEN PG4 = gmael(PG, 1, 4);
    2538           0 :   long lP = lg(PG1);
    2539           0 :   GEN PG5 = pcgrp_lift(gmael(PG, 1, 5), deg);
    2540           0 :   GEN res = cgetg(6, t_VEC), res1, res2, res3;
    2541           0 :   gel(res,1) = res1 = cgetg(lP + 1, t_VEC);
    2542           0 :   gel(res,2) = res2 = cgetg(lP + 1, t_VEC);
    2543           0 :   gel(res,3) = res3 = cgetg(lP + 1, t_VEC);
    2544           0 :   gel(res,4) = vecsmall_prepend(PG4, p);
    2545           0 :   gel(res,5) = PG5;
    2546           0 :   gel(res1, 1) = frob;
    2547           0 :   gel(res2, 1) = ZX_to_Flx(gel(gf->Tmod,1), p);
    2548           0 :   gel(res3, 1) = sigma;
    2549           0 :   for (j = 1; j < lP; j++)
    2550             :   {
    2551             :     struct galois_lift gl;
    2552           0 :     GEN Lden = makeLden(L,den,gb);
    2553             :     GEN pf;
    2554           0 :     initlift(T, den, uel(PG4,j), L, Lden, gb, &gl);
    2555           0 :     pf = nilp_froblift(vecslice(res1,1,j), PG3, j, PG5, gel(PG2,j), incl, H, &gl, gp);
    2556           0 :     if (!pf) return NULL;
    2557           0 :     if (DEBUGLEVEL>=2)
    2558           0 :       err_printf("found: %ld/%ld: %Ps: %Ps\n", n, j+1, gel(pf,2),gel(pf,1));
    2559           0 :     pcgrp_insert(PG5, j+1, gel(pf,2));
    2560           0 :     gel(res1, j+1) = gel(pf,1);
    2561           0 :     gel(res2, j+1) = gel(pf,3);
    2562           0 :     gel(res3, j+1) = gcopy(permtoaut(gel(pf,1), gp));
    2563             :   }
    2564           0 :   if (DEBUGLEVEL >= 4) err_printf("GaloisConj: Fini!\n");
    2565           0 :   return res;
    2566             : }
    2567             : 
    2568             : static GEN
    2569        2142 : galoisgenlift(GEN PG, GEN Pg, GEN O, GEN L, GEN M, GEN frob,
    2570             :               struct galois_borne *gb, struct galois_frobenius *gf)
    2571             : {
    2572             :   struct galois_test td;
    2573             :   GEN res, res1;
    2574        2142 :   GEN PG1 = gel(PG, 1), PG2 = gel(PG, 2);
    2575        2142 :   long lP = lg(PG1), j, n = lg(L)-1;
    2576        2142 :   inittest(L, M, gb->bornesol, gb->ladicsol, &td);
    2577        2142 :   res = cgetg(3, t_VEC);
    2578        2142 :   gel(res,1) = res1 = cgetg(lP + 1, t_VEC);
    2579        2142 :   gel(res,2) = vecsmall_prepend(PG2, gf->deg);
    2580        2142 :   gel(res1, 1) = vecsmall_copy(frob);
    2581        5292 :   for (j = 1; j < lP; j++)
    2582             :   {
    2583        3171 :     GEN pf = galoisgenliftauto(O, gel(PG1, j), gf->psi[Pg[j]], n, &td);
    2584        3171 :     if (!pf) { freetest(&td); return NULL; }
    2585        3150 :     gel(res1, j+1) = pf;
    2586             :   }
    2587        2121 :   if (DEBUGLEVEL >= 4) err_printf("GaloisConj: Fini!\n");
    2588        2121 :   freetest(&td);
    2589        2121 :   return res;
    2590             : }
    2591             : 
    2592             : static ulong
    2593        3024 : psi_order(GEN psi, ulong d)
    2594             : {
    2595        3024 :   long i, l = lg(psi);
    2596        3024 :   ulong s = 1;
    2597        9898 :   for (i=1; i<l; i++)
    2598        6874 :     s = clcm(s, d/cgcd(uel(psi, i)-1, d));
    2599        3024 :   return s;
    2600             : }
    2601             : 
    2602             : static GEN
    2603        3101 : galoisgen(GEN T, GEN L, GEN M, GEN den, GEN bad, struct galois_borne *gb,
    2604             :           const struct galois_analysis *ga)
    2605             : {
    2606             :   struct galois_test td;
    2607             :   struct galois_frobenius gf, ogf;
    2608        3101 :   pari_sp ltop = avma;
    2609        3101 :   long x, n = degpol(T), is_central;
    2610             :   long po;
    2611        3101 :   GEN sigma, res, frob, O, PG, V, ofrob = NULL;
    2612             : 
    2613        3101 :   if (!ga->deg) return NULL;
    2614        3101 :   x = varn(T);
    2615        3101 :   if (DEBUGLEVEL >= 9) err_printf("GaloisConj: denominator:%Ps\n", den);
    2616        3101 :   if (n == 12 && ga->ord==3 && !ga->p4)
    2617             :   { /* A4 is very probable: test it first */
    2618          42 :     pari_sp av = avma;
    2619          42 :     if (DEBUGLEVEL >= 4) err_printf("GaloisConj: Testing A4 first\n");
    2620          42 :     inittest(L, M, gb->bornesol, gb->ladicsol, &td);
    2621          42 :     PG = a4galoisgen(&td);
    2622          42 :     freetest(&td);
    2623          42 :     if (PG) return gerepileupto(ltop, PG);
    2624           0 :     set_avma(av);
    2625             :   }
    2626        3059 :   if (n == 24 && ga->ord==3)
    2627             :   { /* S4 is very probable: test it first */
    2628          35 :     pari_sp av = avma;
    2629             :     struct galois_lift gl;
    2630          35 :     if (DEBUGLEVEL >= 4) err_printf("GaloisConj: Testing S4 first\n");
    2631          35 :     initlift(T, den, ga->p4, L, makeLden(L,den,gb), gb, &gl);
    2632          35 :     PG = s4galoisgen(&gl);
    2633          35 :     if (PG) return gerepileupto(ltop, PG);
    2634           0 :     set_avma(av);
    2635             :   }
    2636        3024 :   frob = galoisfindfrobenius(T, L, den, bad, &gf, gb, ga);
    2637        3024 :   if (!frob) return gc_NULL(ltop);
    2638        3024 :   po = psi_order(gf.psi, gf.deg);
    2639        3024 :   if (!(ga->group&ga_easy) && po < gf.deg && gf.deg/radicalu(gf.deg)%po == 0)
    2640             :   {
    2641           0 :     is_central = 1;
    2642           0 :     if (!bad) bad = gb->dis;
    2643           0 :     if (po > 1)
    2644             :     {
    2645           0 :       ofrob = frob; ogf = gf;
    2646           0 :       frob = perm_pow(frob, po);
    2647           0 :       gf.deg /= po;
    2648             :     }
    2649        3024 :   } else is_central = 0;
    2650        3024 :   sigma = permtopol(frob, L, M, den, gb->ladicabs, shifti(gb->ladicabs,-1), x);
    2651        3024 :   if (is_central && gf.fp != gf.deg)
    2652           0 :   { gf.p = findp(bad, T, sigma, gf.deg, &gf.Tmod); gf.fp = gf.deg;
    2653           0 :     gf.psi = const_vecsmall(lg(gf.Tmod)-1, 1);
    2654             :   }
    2655        3024 :   if (gf.deg == n)        /* cyclic */
    2656             :   {
    2657         861 :     GEN Tp = ZX_to_Flx(gel(gf.Tmod,1), gf.p);
    2658         861 :     res = mkvec5(mkvec(frob), mkvec(Tp), mkvec(sigma), mkvecsmall(gf.p), cyclic_pc(n));
    2659         861 :     return gerepilecopy(ltop, res);
    2660             :   }
    2661        2163 :   O = perm_cycles(frob);
    2662        2163 :   if (DEBUGLEVEL >= 9) err_printf("GaloisConj: Frobenius:%Ps\n", sigma);
    2663        2163 :   PG = galoisgenfixedfield0(O, L, sigma, T, is_central ? bad: NULL,
    2664             :                                             is_central ? &V:  NULL, &gf, gb);
    2665        2163 :   if (PG == NULL) return gc_NULL(ltop);
    2666        2142 :   if (is_central && lg(gel(PG,1))!=3)
    2667           0 :   {
    2668             :     struct galois_perm gp;
    2669           0 :     galoisperm_init(&gp, L, M, den, gb->ladicabs, shifti(gb->ladicabs,-1), varn(T));
    2670           0 :     res = galoisgenlift_nilp(PG, O, V, T, frob, sigma, gb, &gf, &gp);
    2671           0 :     galoisperm_free(&gp);
    2672             :   }
    2673             :   else
    2674             :   {
    2675        2142 :     if (is_central && po > 1)
    2676             :     { /* backtrack powering of frob */
    2677           0 :       frob = ofrob; gf = ogf;
    2678           0 :       O = perm_cycles(ofrob);
    2679           0 :       sigma = permtopol(ofrob, L, M, den, gb->ladicabs, shifti(gb->ladicabs,-1), x);
    2680           0 :       PG = galoisgenfixedfield0(O, L, sigma, T, NULL, NULL, &gf, gb);
    2681           0 :       if (PG == NULL) return gc_NULL(ltop);
    2682             :     }
    2683        2142 :     res = galoisgenlift(gg_get_std(gel(PG,1)), gel(PG,2), O, L, M, frob, gb, &gf);
    2684             :   }
    2685        2142 :   if (!res) return gc_NULL(ltop);
    2686        2121 :   return gerepilecopy(ltop, res);
    2687             : }
    2688             : 
    2689             : /* T = polcyclo(N) */
    2690             : static GEN
    2691          70 : conjcyclo(GEN T, long N)
    2692             : {
    2693          70 :   pari_sp av = avma;
    2694          70 :   long i, k = 1, d = eulerphiu(N), v = varn(T);
    2695          70 :   GEN L = cgetg(d+1,t_COL);
    2696         651 :   for (i=1; i<=N; i++)
    2697         581 :     if (ugcd(i, N)==1)
    2698             :     {
    2699         294 :       GEN s = pol_xn(i, v);
    2700         294 :       if (i >= d) s = ZX_rem(s, T);
    2701         294 :       gel(L,k++) = s;
    2702             :     }
    2703          70 :   return gerepileupto(av, gen_sort(L, (void*)&gcmp, &gen_cmp_RgX));
    2704             : }
    2705             : 
    2706             : static GEN
    2707           0 : aut_to_groupelts(GEN aut, GEN L, ulong p)
    2708             : {
    2709           0 :   pari_sp av = avma;
    2710           0 :   long i, d = lg(aut)-1;
    2711           0 :   GEN P = ZV_to_Flv(L, p);
    2712           0 :   GEN N = FlxV_Flv_multieval(RgXV_to_FlxV(aut, p), P, p);
    2713           0 :   GEN q = perm_inv(vecsmall_indexsort(P));
    2714           0 :   GEN G = cgetg(d+1, t_VEC);
    2715           0 :   for (i=1; i<=d; i++)
    2716           0 :     gel(G,i) = perm_mul(vecsmall_indexsort(gel(N,i)),q);
    2717           0 :   return gerepilecopy(av, vecvecsmall_sort(G));
    2718             : }
    2719             : 
    2720             : static GEN
    2721           0 : galoisinitfromaut(GEN T, GEN aut)
    2722             : {
    2723           0 :   pari_sp ltop = avma;
    2724           0 :   GEN nf, A, G, L, M, grp, den=NULL;
    2725             :   struct galois_analysis ga;
    2726             :   struct galois_borne gb;
    2727             :   long n;
    2728             :   pari_timer ti;
    2729             : 
    2730           0 :   T = get_nfpol(T, &nf);
    2731           0 :   n = degpol(T);
    2732           0 :   if (nf)
    2733           0 :   { if (!den) den = nf_get_zkden(nf); }
    2734             :   else
    2735             :   {
    2736           0 :     if (n <= 0) pari_err_IRREDPOL("galoisinit",T);
    2737           0 :     RgX_check_ZX(T, "galoisinit");
    2738           0 :     if (!ZX_is_squarefree(T))
    2739           0 :       pari_err_DOMAIN("galoisinit","issquarefree(pol)","=",gen_0,T);
    2740           0 :     if (!gequal1(gel(T,n+2))) pari_err_IMPL("galoisinit(non-monic)");
    2741             :   }
    2742           0 :   if (!galoisanalysis(T, &ga, 1, NULL)) pari_err_IMPL("galoisinit");
    2743           0 :   gb.l = utoipos(ga.l);
    2744           0 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
    2745           0 :   den = galoisborne(T, den, &gb, degpol(T));
    2746           0 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "galoisborne()");
    2747           0 :   L = ZpX_roots(T, gb.l, gb.valabs);
    2748           0 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "ZpX_roots");
    2749           0 :   M = FpV_invVandermonde(L, den, gb.ladicabs);
    2750           0 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "FpV_invVandermonde()");
    2751           0 :   A = aut_to_groupelts(aut, L, ga.l);
    2752           0 :   G = groupelts_to_group(A);
    2753           0 :   if (!G) pari_err_IMPL("galoisinit");
    2754           0 :   A = group_elts(G,n);
    2755           0 :   grp = cgetg(9, t_VEC);
    2756           0 :   gel(grp,1) = T;
    2757           0 :   gel(grp,2) = mkvec3(utoipos(ga.l), utoipos(gb.valabs), gb.ladicabs);
    2758           0 :   gel(grp,3) = L;
    2759           0 :   gel(grp,4) = M;
    2760           0 :   gel(grp,5) = den;
    2761           0 :   gel(grp,6) = A;
    2762           0 :   gel(grp,7) = gel(G,1);
    2763           0 :   gel(grp,8) = gel(G,2);
    2764           0 :   return gerepilecopy(ltop, grp);
    2765             : }
    2766             : 
    2767             : /* T: polynomial or nf, den multiple of common denominator of solutions or
    2768             :  * NULL (unknown). If T is nf, and den unknown, use den = denom(nf.zk) */
    2769             : static GEN
    2770        4501 : galoisconj4_main(GEN T, GEN den, long flag)
    2771             : {
    2772        4501 :   pari_sp ltop = avma;
    2773             :   GEN nf, G, L, M, aut, grp;
    2774             :   struct galois_analysis ga;
    2775             :   struct galois_borne gb;
    2776             :   long n;
    2777             :   pari_timer ti;
    2778             : 
    2779        4501 :   T = get_nfpol(T, &nf);
    2780        4501 :   n = poliscyclo(T);
    2781        4501 :   if (n) return flag? galoiscyclo(n, varn(T)): conjcyclo(T, n);
    2782        4207 :   n = degpol(T);
    2783        4207 :   if (nf)
    2784        3171 :   { if (!den) den = nf_get_zkden(nf); }
    2785             :   else
    2786             :   {
    2787        1036 :     if (n <= 0) pari_err_IRREDPOL("galoisinit",T);
    2788        1036 :     RgX_check_ZX(T, "galoisinit");
    2789        1036 :     if (!ZX_is_squarefree(T))
    2790           7 :       pari_err_DOMAIN("galoisinit","issquarefree(pol)","=",gen_0,T);
    2791        1029 :     if (!gequal1(gel(T,n+2))) pari_err_IMPL("galoisinit(non-monic)");
    2792             :   }
    2793        4193 :   if (n == 1)
    2794             :   {
    2795          21 :     if (!flag) { G = cgetg(2, t_COL); gel(G,1) = pol_x(varn(T)); return G;}
    2796          21 :     ga.l = 3;
    2797          21 :     ga.deg = 1;
    2798          21 :     den = gen_1;
    2799             :   }
    2800        4172 :   else if (!galoisanalysis(T, &ga, 1, NULL)) return gc_NULL(ltop);
    2801             : 
    2802        2366 :   if (den)
    2803             :   {
    2804        1841 :     if (typ(den) != t_INT) pari_err_TYPE("galoisconj4 [2nd arg integer]", den);
    2805        1841 :     den = absi_shallow(den);
    2806             :   }
    2807        2366 :   gb.l = utoipos(ga.l);
    2808        2366 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
    2809        2366 :   den = galoisborne(T, den, &gb, degpol(T));
    2810        2366 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "galoisborne()");
    2811        2366 :   L = ZpX_roots(T, gb.l, gb.valabs);
    2812        2366 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "ZpX_roots");
    2813        2366 :   M = FpV_invVandermonde(L, den, gb.ladicabs);
    2814        2366 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "FpV_invVandermonde()");
    2815        2366 :   if (n == 1)
    2816             :   {
    2817          21 :     G = cgetg(3, t_VEC);
    2818          21 :     gel(G,1) = cgetg(1, t_VEC);
    2819          21 :     gel(G,2) = cgetg(1, t_VECSMALL);
    2820             :   }
    2821             :   else
    2822        2345 :     G = gg_get_std(galoisgen(T, L, M, den, NULL, &gb, &ga));
    2823        2366 :   if (DEBUGLEVEL >= 6) err_printf("GaloisConj: %Ps\n", G);
    2824        2366 :   if (!G) return gc_NULL(ltop);
    2825        2331 :   if (DEBUGLEVEL >= 1) timer_start(&ti);
    2826        2331 :   grp = cgetg(9, t_VEC);
    2827        2331 :   gel(grp,1) = T;
    2828        2331 :   gel(grp,2) = mkvec3(utoipos(ga.l), utoipos(gb.valabs), gb.ladicabs);
    2829        2331 :   gel(grp,3) = L;
    2830        2331 :   gel(grp,4) = M;
    2831        2331 :   gel(grp,5) = den;
    2832        2331 :   gel(grp,6) = group_elts(G,n);
    2833        2331 :   gel(grp,7) = gel(G,1);
    2834        2331 :   gel(grp,8) = gel(G,2);
    2835        2331 :   if (flag) return gerepilecopy(ltop, grp);
    2836         945 :   aut = galoisvecpermtopol(grp, gal_get_group(grp), gb.ladicabs, shifti(gb.ladicabs,-1));
    2837         945 :   settyp(aut, t_COL);
    2838         945 :   if (DEBUGLEVEL >= 1) timer_printf(&ti, "Computation of polynomials");
    2839         945 :   return gerepileupto(ltop, gen_sort(aut, (void*)&gcmp, &gen_cmp_RgX));
    2840             : }
    2841             : 
    2842             : /* Heuristic computation of #Aut(T), pinit = first prime to be tested */
    2843             : long
    2844         917 : numberofconjugates(GEN T, long pinit)
    2845             : {
    2846         917 :   pari_sp av = avma;
    2847         917 :   long c, nbtest, nbmax, n = degpol(T);
    2848             :   ulong p;
    2849             :   forprime_t S;
    2850             : 
    2851         917 :   if (n == 1) return 1;
    2852         917 :   nbmax = (n < 10)? 20: (n<<1) + 1;
    2853         917 :   nbtest = 0;
    2854             : #if 0
    2855             :   c = ZX_sturm(T); c = ugcd(c, n-c); /* too costly: finite primes are cheaper */
    2856             : #else
    2857         917 :   c = n;
    2858             : #endif
    2859         917 :   u_forprime_init(&S, pinit, ULONG_MAX);
    2860       10472 :   while((p = u_forprime_next(&S)))
    2861             :   {
    2862       10472 :     GEN L, Tp = ZX_to_Flx(T,p);
    2863             :     long i, nb;
    2864       10472 :     if (!Flx_is_squarefree(Tp, p)) continue;
    2865             :     /* unramified */
    2866        8841 :     nbtest++;
    2867        8841 :     L = Flx_nbfact_by_degree(Tp, &nb, p); /* L[i] = #factors of degree i */
    2868        8841 :     if (L[n/nb] == nb) {
    2869        6503 :       if (c == n && nbtest > 10) break; /* probably Galois */
    2870             :     }
    2871             :     else
    2872             :     {
    2873        3220 :       c = ugcd(c, L[1]);
    2874       30058 :       for (i = 2; i <= n; i++)
    2875       27377 :         if (L[i]) { c = ugcd(c, L[i]*i); if (c == 1) break; }
    2876        3220 :       if (c == 1) break;
    2877             :     }
    2878        8267 :     if (nbtest == nbmax) break;
    2879        7924 :     if (DEBUGLEVEL >= 6)
    2880           0 :       err_printf("NumberOfConjugates [%ld]:c=%ld,p=%ld\n", nbtest,c,p);
    2881        7924 :     set_avma(av);
    2882             :   }
    2883         917 :   if (DEBUGLEVEL >= 2) err_printf("NumberOfConjugates:c=%ld,p=%ld\n", c, p);
    2884         917 :   return gc_long(av,c);
    2885             : }
    2886             : static GEN
    2887           0 : galoisconj4(GEN nf, GEN d)
    2888             : {
    2889           0 :   pari_sp av = avma;
    2890             :   GEN G, T;
    2891           0 :   G = galoisconj4_main(nf, d, 0);
    2892           0 :   if (G) return G; /* Success */
    2893           0 :   set_avma(av); T = get_nfpol(nf, &nf);
    2894           0 :   G = cgetg(2, t_COL); gel(G,1) = pol_x(varn(T)); return G; /* Fail */
    2895             : 
    2896             : }
    2897             : 
    2898             : /* d multiplicative bound for the automorphism's denominators */
    2899             : GEN
    2900       12257 : galoisconj(GEN nf, GEN d)
    2901             : {
    2902       12257 :   pari_sp av = avma;
    2903       12257 :   GEN G, NF, T = get_nfpol(nf,&NF);
    2904       12257 :   if (degpol(T) == 2)
    2905             :   { /* fast shortcut */
    2906       10332 :     GEN a = gel(T,4), b = gel(T,3);
    2907       10332 :     long v = varn(T);
    2908       10332 :     RgX_check_ZX(T, "galoisconj");
    2909       10332 :     if (!gequal1(a)) pari_err_IMPL("galoisconj(non-monic)");
    2910       10332 :     G = cgetg(3, t_COL);
    2911       10332 :     gel(G,1) = deg1pol_shallow(gen_m1, negi(b), v);
    2912       10332 :     gel(G,2) = pol_x(v);
    2913       10332 :     return G;
    2914             :   }
    2915        1925 :   G = galoisconj4_main(nf, d, 0);
    2916        1925 :   if (G) return G; /* Success */
    2917         910 :   set_avma(av); return galoisconj1(nf);
    2918             : }
    2919             : 
    2920             : /* FIXME: obsolete, use galoisconj(nf, d) directly */
    2921             : GEN
    2922          49 : galoisconj0(GEN nf, long flag, GEN d, long prec)
    2923             : {
    2924             :   (void)prec;
    2925          49 :   switch(flag) {
    2926          42 :     case 2:
    2927          42 :     case 0: return galoisconj(nf, d);
    2928           7 :     case 1: return galoisconj1(nf);
    2929           0 :     case 4: return galoisconj4(nf, d);
    2930             :   }
    2931           0 :   pari_err_FLAG("nfgaloisconj");
    2932             :   return NULL; /*LCOV_EXCL_LINE*/
    2933             : }
    2934             : 
    2935             : /******************************************************************************/
    2936             : /* Galois theory related algorithms                                           */
    2937             : /******************************************************************************/
    2938             : GEN
    2939       20699 : checkgal(GEN gal)
    2940             : {
    2941       20699 :   if (typ(gal) == t_POL) pari_err_TYPE("checkgal [apply galoisinit first]",gal);
    2942       20699 :   if (typ(gal) != t_VEC || lg(gal) != 9) pari_err_TYPE("checkgal",gal);
    2943       20692 :   return gal;
    2944             : }
    2945             : 
    2946             : GEN
    2947        2576 : galoisinit(GEN nf, GEN den)
    2948             : {
    2949             :   GEN G;
    2950        2576 :   if (is_vec_t(typ(nf)) && lg(nf)==3 && is_vec_t(typ(gel(nf,2))))
    2951           0 :     return galoisinitfromaut(gel(nf,1), gel(nf,2));
    2952        2576 :   G = galoisconj4_main(nf, den, 1);
    2953        2562 :   return G? G: gen_0;
    2954             : }
    2955             : 
    2956             : static GEN
    2957       14952 : galoispermtopol_i(GEN gal, GEN perm, GEN mod, GEN mod2)
    2958             : {
    2959       14952 :   switch (typ(perm))
    2960             :   {
    2961       14714 :     case t_VECSMALL:
    2962       14714 :       return permtopol(perm, gal_get_roots(gal), gal_get_invvdm(gal),
    2963             :                              gal_get_den(gal), mod, mod2,
    2964       14714 :                              varn(gal_get_pol(gal)));
    2965         238 :     case t_VEC: case t_COL: case t_MAT:
    2966         238 :       return galoisvecpermtopol(gal, perm, mod, mod2);
    2967             :   }
    2968           0 :   pari_err_TYPE("galoispermtopol", perm);
    2969             :   return NULL; /* LCOV_EXCL_LINE */
    2970             : }
    2971             : 
    2972             : GEN
    2973       14952 : galoispermtopol(GEN gal, GEN perm)
    2974             : {
    2975       14952 :   pari_sp av = avma;
    2976             :   GEN mod, mod2;
    2977       14952 :   gal = checkgal(gal);
    2978       14952 :   mod = gal_get_mod(gal);
    2979       14952 :   mod2 = shifti(mod,-1);
    2980       14952 :   return gerepilecopy(av, galoispermtopol_i(gal, perm, mod, mod2));
    2981             : }
    2982             : 
    2983             : GEN
    2984         119 : galoiscosets(GEN O, GEN perm)
    2985             : {
    2986         119 :   long i, j, k, u, f, l = lg(O);
    2987         119 :   GEN RC, C = cgetg(l,t_VECSMALL), o = gel(O,1);
    2988         119 :   pari_sp av = avma;
    2989         119 :   f = lg(o); u = o[1]; RC = zero_zv(lg(perm)-1);
    2990         546 :   for(i=1,j=1; j<l; i++)
    2991             :   {
    2992         427 :     GEN p = gel(perm,i);
    2993         427 :     if (RC[ p[u] ]) continue;
    2994        1001 :     for(k=1; k<f; k++) RC[ p[ o[k] ] ] = 1;
    2995         357 :     C[j++] = i;
    2996             :   }
    2997         119 :   set_avma(av); return C;
    2998             : }
    2999             : 
    3000             : static GEN
    3001         119 : fixedfieldfactor(GEN L, GEN O, GEN perm, GEN M, GEN den, GEN mod, GEN mod2,
    3002             :                  long x,long y)
    3003             : {
    3004         119 :   pari_sp ltop = avma;
    3005         119 :   long i, j, k, l = lg(O), lo = lg(gel(O,1));
    3006         119 :   GEN V, res, cosets = galoiscosets(O,perm), F = cgetg(lo+1,t_COL);
    3007             : 
    3008         119 :   gel(F, lo) = gen_1;
    3009         119 :   if (DEBUGLEVEL>=4) err_printf("GaloisFixedField:cosets=%Ps \n",cosets);
    3010         119 :   if (DEBUGLEVEL>=6) err_printf("GaloisFixedField:den=%Ps mod=%Ps \n",den,mod);
    3011         119 :   V = cgetg(l,t_COL); res = cgetg(l,t_VEC);
    3012         476 :   for (i = 1; i < l; i++)
    3013             :   {
    3014         357 :     pari_sp av = avma;
    3015         357 :     GEN G = cgetg(l,t_VEC), Lp = vecpermute(L, gel(perm, cosets[i]));
    3016        1680 :     for (k = 1; k < l; k++)
    3017        1323 :       gel(G,k) = FpV_roots_to_pol(vecpermute(Lp, gel(O,k)), mod, x);
    3018        1001 :     for (j = 1; j < lo; j++)
    3019             :     {
    3020        2548 :       for(k = 1; k < l; k++) gel(V,k) = gmael(G,k,j+1);
    3021         644 :       gel(F,j) = vectopol(V, M, den, mod, mod2, y);
    3022             :     }
    3023         357 :     gel(res,i) = gerepileupto(av,gtopolyrev(F,x));
    3024             :   }
    3025         119 :   return gerepileupto(ltop,res);
    3026             : }
    3027             : 
    3028             : static void
    3029        1701 : chk_perm(GEN perm, long n)
    3030             : {
    3031        1701 :   if (typ(perm) != t_VECSMALL || lg(perm)!=n+1)
    3032           0 :     pari_err_TYPE("galoisfixedfield", perm);
    3033        1701 : }
    3034             : 
    3035             : static int
    3036        6279 : is_group(GEN g)
    3037             : {
    3038        6272 :   return typ(g)==t_VEC && lg(g)==3 && typ(gel(g,1))==t_VEC
    3039       12551 :       && typ(gel(g,2))==t_VECSMALL;
    3040             : }
    3041             : 
    3042             : GEN
    3043        1211 : galoisfixedfield(GEN gal, GEN perm, long flag, long y)
    3044             : {
    3045        1211 :   pari_sp ltop = avma;
    3046             :   GEN T, L, P, S, PL, O, res, mod, mod2, OL, sym;
    3047             :   long vT, n, i;
    3048        1211 :   if (flag<0 || flag>2) pari_err_FLAG("galoisfixedfield");
    3049        1211 :   gal = checkgal(gal); T = gal_get_pol(gal);
    3050        1211 :   vT = varn(T);
    3051        1211 :   L = gal_get_roots(gal); n = lg(L)-1;
    3052        1211 :   mod = gal_get_mod(gal);
    3053        1211 :   if (typ(perm) == t_VEC)
    3054             :   {
    3055         854 :     if (is_group(perm)) perm = gel(perm, 1);
    3056        2198 :     for (i = 1; i < lg(perm); i++) chk_perm(gel(perm,i), n);
    3057         854 :     O = vecperm_orbits(perm, n);
    3058             :   }
    3059             :   else
    3060             :   {
    3061         357 :     chk_perm(perm, n);
    3062         357 :     O = perm_cycles(perm);
    3063             :   }
    3064        1211 :   mod2 = shifti(mod,-1);
    3065        1211 :   OL = fixedfieldorbits(O, L);
    3066        1211 :   sym = fixedfieldsympol(OL, itou(gal_get_p(gal)));
    3067        1211 :   PL = sympol_eval(sym, OL, mod);
    3068        1211 :   P = FpX_center_i(FpV_roots_to_pol(PL, mod, vT), mod, mod2);
    3069        1211 :   if (flag==1) return gerepilecopy(ltop,P);
    3070         875 :   S = fixedfieldinclusion(O, PL);
    3071         875 :   S = vectopol(S, gal_get_invvdm(gal), gal_get_den(gal), mod, mod2, vT);
    3072         875 :   if (flag==0)
    3073         756 :     res = cgetg(3, t_VEC);
    3074             :   else
    3075             :   {
    3076             :     GEN PM, Pden;
    3077             :     struct galois_borne Pgb;
    3078         119 :     long val = itos(gal_get_e(gal));
    3079         119 :     Pgb.l = gal_get_p(gal);
    3080         119 :     Pden = galoisborne(P, NULL, &Pgb, degpol(T)/degpol(P));
    3081         119 :     if (Pgb.valabs > val)
    3082             :     {
    3083          20 :       if (DEBUGLEVEL>=4)
    3084           0 :         err_printf("GaloisConj: increase p-adic prec by %ld.\n", Pgb.valabs-val);
    3085          20 :       PL = ZpX_liftroots(P, PL, Pgb.l, Pgb.valabs);
    3086          20 :       L  = ZpX_liftroots(T, L, Pgb.l, Pgb.valabs);
    3087          20 :       mod = Pgb.ladicabs; mod2 = shifti(mod,-1);
    3088             :     }
    3089         119 :     PM = FpV_invVandermonde(PL, Pden, mod);
    3090         119 :     if (y < 0) y = 1;
    3091         119 :     if (varncmp(y, vT) <= 0)
    3092           0 :       pari_err_PRIORITY("galoisfixedfield", T, "<=", y);
    3093         119 :     setvarn(P, y);
    3094         119 :     res = cgetg(4, t_VEC);
    3095         119 :     gel(res,3) = fixedfieldfactor(L,O,gal_get_group(gal), PM,Pden,mod,mod2,vT,y);
    3096             :   }
    3097         875 :   gel(res,1) = gcopy(P);
    3098         875 :   gel(res,2) = gmodulo(S, T);
    3099         875 :   return gerepileupto(ltop, res);
    3100             : }
    3101             : 
    3102             : /* gal a galois group output the underlying wss group */
    3103             : GEN
    3104        1330 : galois_group(GEN gal) { return mkvec2(gal_get_gen(gal), gal_get_orders(gal)); }
    3105             : 
    3106             : GEN
    3107         826 : checkgroup(GEN g, GEN *S)
    3108             : {
    3109         826 :   if (is_group(g)) { *S = NULL; return g; }
    3110         469 :   g  = checkgal(g);
    3111         462 :   *S = gal_get_group(g); return galois_group(g);
    3112             : }
    3113             : 
    3114             : GEN
    3115        4613 : checkgroupelts(GEN G)
    3116             : {
    3117             :   long i, n;
    3118        4613 :   if (typ(G)!=t_VEC) pari_err_TYPE("checkgroupelts", G);
    3119        4599 :   if (is_group(G))
    3120             :   { /* subgroup of S_n */
    3121         371 :     if (lg(gel(G,1))==1) return mkvec(mkvecsmall(1));
    3122         371 :     return group_elts(G, group_domain(G));
    3123             :   }
    3124        4228 :   if (lg(G)==9 && typ(gel(G,1))==t_POL)
    3125        3913 :     return gal_get_group(G); /* galoisinit */
    3126             :   /* vector of permutations ? */
    3127         315 :   n = lg(G)-1;
    3128         315 :   if (n==0) pari_err_DIM("checkgroupelts");
    3129        4984 :   for (i = 1; i <= n; i++)
    3130             :   {
    3131        4704 :     if (typ(gel(G,i)) != t_VECSMALL)
    3132          21 :       pari_err_TYPE("checkgroupelts (element)", gel(G,i));
    3133        4683 :     if (lg(gel(G,i)) != lg(gel(G,1)))
    3134          14 :       pari_err_DIM("checkgroupelts [length of permutations]");
    3135             :   }
    3136         280 :   return G;
    3137             : }
    3138             : 
    3139             : GEN
    3140         168 : galoisisabelian(GEN gal, long flag)
    3141             : {
    3142         168 :   pari_sp av = avma;
    3143         168 :   GEN S, G = checkgroup(gal,&S);
    3144         168 :   if (!group_isabelian(G)) { set_avma(av); return gen_0; }
    3145         147 :   switch(flag)
    3146             :   {
    3147          49 :     case 0: return gerepileupto(av, group_abelianHNF(G,S));
    3148          49 :     case 1: set_avma(av); return gen_1;
    3149          49 :     case 2: return gerepileupto(av, group_abelianSNF(G,S));
    3150           0 :     default: pari_err_FLAG("galoisisabelian");
    3151             :   }
    3152             :   return NULL; /* LCOV_EXCL_LINE */
    3153             : }
    3154             : 
    3155             : long
    3156          56 : galoisisnormal(GEN gal, GEN sub)
    3157             : {
    3158          56 :   pari_sp av = avma;
    3159          56 :   GEN S, G = checkgroup(gal, &S), H = checkgroup(sub, &S);
    3160          56 :   long res = group_subgroup_isnormal(G, H);
    3161          56 :   set_avma(av);
    3162          56 :   return res;
    3163             : }
    3164             : 
    3165             : static GEN
    3166         308 : conjclasses_count(GEN conj, long nb)
    3167             : {
    3168         308 :   long i, l = lg(conj);
    3169         308 :   GEN c = zero_zv(nb);
    3170        4039 :   for (i = 1; i < l; i++) c[conj[i]]++;
    3171         308 :   return c;
    3172             : }
    3173             : GEN
    3174         308 : galoisconjclasses(GEN G)
    3175             : {
    3176         308 :   pari_sp av = avma;
    3177         308 :   GEN c, e, cc = group_to_cc(G);
    3178         308 :   GEN elts = gel(cc,1), conj = gel(cc,2), repr = gel(cc,3);
    3179         308 :   long i, l = lg(conj), lc = lg(repr);
    3180         308 :   c = conjclasses_count(conj, lc-1);
    3181         308 :   e = cgetg(lc, t_VEC);
    3182        3143 :   for (i = 1; i < lc; i++) gel(e,i) = cgetg(c[i]+1, t_VEC);
    3183        4039 :   for (i = 1; i < l; i++)
    3184             :   {
    3185        3731 :     long ci = conj[i];
    3186        3731 :     gmael(e, ci, c[ci]) = gel(elts, i);
    3187        3731 :     c[ci]--;
    3188             :   }
    3189         308 :   return gerepilecopy(av, e);
    3190             : }
    3191             : 
    3192             : GEN
    3193          77 : galoissubgroups(GEN gal)
    3194             : {
    3195          77 :   pari_sp av = avma;
    3196          77 :   GEN S, G = checkgroup(gal,&S);
    3197          77 :   return gerepileupto(av, group_subgroups(G));
    3198             : }
    3199             : 
    3200             : GEN
    3201          56 : galoissubfields(GEN G, long flag, long v)
    3202             : {
    3203          56 :   pari_sp av = avma;
    3204          56 :   GEN L = galoissubgroups(G);
    3205          56 :   long i, l = lg(L);
    3206          56 :   GEN S = cgetg(l, t_VEC);
    3207         847 :   for (i = 1; i < l; ++i) gel(S,i) = galoisfixedfield(G, gmael(L,i,1), flag, v);
    3208          56 :   return gerepileupto(av, S);
    3209             : }
    3210             : 
    3211             : GEN
    3212          28 : galoisexport(GEN gal, long format)
    3213             : {
    3214          28 :   pari_sp av = avma;
    3215          28 :   GEN S, G = checkgroup(gal,&S);
    3216          28 :   return gerepileupto(av, group_export(G,format));
    3217             : }
    3218             : 
    3219             : GEN
    3220         406 : galoisidentify(GEN gal)
    3221             : {
    3222         406 :   pari_sp av = avma;
    3223         406 :   GEN S, G = checkgroup(gal,&S);
    3224         399 :   long idx = group_ident(G,S), card = group_order(G);
    3225         399 :   set_avma(av); return mkvec2s(card, idx);
    3226             : }
    3227             : 
    3228             : /* index of conjugacy class containing g */
    3229             : static long
    3230       36939 : cc_id(GEN cc, GEN g)
    3231             : {
    3232       36939 :   GEN conj = gel(cc,2);
    3233       36939 :   long k = signe(gel(cc,4))? g[1]: vecvecsmall_search(gel(cc,1),g,0);
    3234       36939 :   return conj[k];
    3235             : }
    3236             : 
    3237             : static GEN
    3238        4186 : Qevproj_RgX(GEN c, long d, GEN pro)
    3239        4186 : { return RgV_to_RgX(Qevproj_down(RgX_to_RgC(c,d), pro), varn(c)); }
    3240             : /* c in Z[X] / (X^o-1), To = polcyclo(o), T = polcyclo(expo), e = expo/o
    3241             :  * return c(X^e) mod T as an element of Z[X] / (To) */
    3242             : static GEN
    3243        3920 : chival(GEN c, GEN T, GEN To, long e, GEN pro, long phie)
    3244             : {
    3245        3920 :   c = ZX_rem(c, To);
    3246        3920 :   if (e != 1) c = ZX_rem(RgX_inflate(c,e), T);
    3247        3920 :   if (pro) c = Qevproj_RgX(c, phie, pro);
    3248        3920 :   return c;
    3249             : }
    3250             : /* chi(g^l) = sum_{k=0}^{o-1} a_k zeta_o^{l*k} for all l;
    3251             : * => a_k = 1/o sum_{l=0}^{o-1} chi(g^l) zeta_o^{-k*l}. Assume o > 1 */
    3252             : static GEN
    3253         861 : chiFT(GEN cp, GEN jg, GEN vze, long e, long o, ulong p, ulong pov2)
    3254             : {
    3255         861 :   const long var = 1;
    3256         861 :   ulong oinv = Fl_inv(o,p);
    3257             :   long k, l;
    3258         861 :   GEN c = cgetg(o+2, t_POL);
    3259        5642 :   for (k = 0; k < o; k++)
    3260             :   {
    3261        4781 :     ulong a = 0;
    3262       51478 :     for (l=0; l<o; l++)
    3263             :     {
    3264       46697 :       ulong z = vze[Fl_mul(k,l,o)*e + 1];/* zeta_o^{-k*l} */
    3265       46697 :       a = Fl_add(a, Fl_mul(uel(cp,jg[l+1]), z, p), p);
    3266             :     }
    3267        4781 :     gel(c,k+2) = stoi(Fl_center(Fl_mul(a,oinv,p), p, pov2)); /* a_k */
    3268             :   }
    3269         861 :   c[1] = evalvarn(var) | evalsigne(1); return ZX_renormalize(c,o+2);
    3270             : }
    3271             : static GEN
    3272         546 : cc_chartable(GEN cc)
    3273             : {
    3274             :   GEN al, elts, rep, ctp, ct, dec, id, vjg, H, vord, operm;
    3275             :   long i, j, k, f, l, expo, lcl, n;
    3276             :   ulong p, pov2;
    3277             : 
    3278         546 :   elts = gel(cc,1); n = lg(elts)-1;
    3279         546 :   if (n == 1) return mkvec2(mkmat(mkcol(gen_1)), gen_1);
    3280         532 :   rep = gel(cc,3);
    3281         532 :   lcl = lg(rep);
    3282         532 :   vjg = cgetg(lcl, t_VEC);
    3283         532 :   vord = cgetg(lcl,t_VECSMALL);
    3284         532 :   id = identity_perm(lg(gel(elts,1))-1);
    3285         532 :   expo = 1;
    3286        4879 :   for(j=1;j<lcl;j++)
    3287             :   {
    3288        4347 :     GEN jg, h = id, g = gel(elts,rep[j]);
    3289             :     long o;
    3290        4347 :     vord[j] = o = perm_order(g);
    3291        4347 :     expo = ulcm(expo, o);
    3292        4347 :     gel(vjg,j) = jg = cgetg(o+1,t_VECSMALL);
    3293       27671 :     for (l=1; l<=o; l++)
    3294             :     {
    3295       23324 :       jg[l] = cc_id(cc, h); /* index of conjugacy class of g^(l-1) */
    3296       23324 :       if (l < o) h = perm_mul(h, g);
    3297             :     }
    3298             :   }
    3299             :   /* would sort conjugacy classes by inc. order */
    3300         532 :   operm = vecsmall_indexsort(vord);
    3301             : 
    3302             :   /* expo > 1, exponent of G */
    3303         532 :   p = unextprime(2*n+1);
    3304        1043 :   while (p%expo != 1) p = unextprime(p+1);
    3305             :   /* compute character table modulo p: idempotents of Z(KG) */
    3306         532 :   al = conjclasses_algcenter(cc, utoipos(p));
    3307         532 :   dec = algsimpledec_ss(al,1);
    3308         532 :   ctp = cgetg(lcl,t_VEC);
    3309        4879 :   for(i=1; i<lcl; i++)
    3310             :   {
    3311        4347 :     GEN e = ZV_to_Flv(gmael3(dec,i,3,1), p); /*(1/n)[(dim chi)chi(g): g in G]*/
    3312        4347 :     ulong d = usqrt(Fl_mul(e[1], n, p)); /* = chi(1) <= sqrt(n) < sqrt(p) */
    3313        4347 :     gel(ctp,i) = Flv_Fl_mul(e,Fl_div(n,d,p), p); /*[chi(g): g in G]*/
    3314             :   }
    3315             :   /* Find minimal f such that table is defined over Q(zeta(f)): the conductor
    3316             :    * of the class field Q(\zeta_e)^H defined by subgroup
    3317             :    * H = { k in (Z/e)^*: g^k ~ g, for all g } */
    3318         532 :   H = coprimes_zv(expo);
    3319        3458 :   for (k = 2; k < expo; k++)
    3320             :   {
    3321        2926 :     if (!H[k]) continue;
    3322        2548 :     for (j = 2; j < lcl; j++) /* skip g ~ 1 */
    3323        2366 :       if (umael(vjg,j,(k % vord[j])+1) != umael(vjg,j,2)) { H[k] = 0; break; }
    3324             :   }
    3325         532 :   f = znstar_conductor_bits(Flv_to_F2v(H));
    3326             :   /* lift character table to Z[zeta_f] */
    3327         532 :   pov2 = p>>1;
    3328         532 :   ct = cgetg(lcl, t_MAT);
    3329         532 :   if (f == 1)
    3330             :   { /* rational representation */
    3331         938 :     for (j=1; j<lcl; j++) gel(ct,j) = cgetg(lcl,t_COL);
    3332         938 :     for(j=1; j<lcl; j++)
    3333             :     {
    3334         791 :       GEN jg = gel(vjg,j); /* jg[l+1] = class of g^l */
    3335         791 :       long t = lg(jg) > 2? jg[2]: jg[1];
    3336        6706 :       for(i=1; i<lcl; i++)
    3337             :       {
    3338        5915 :         GEN cp = gel(ctp,i); /* cp[i] = chi(g_i) mod \P */
    3339        5915 :         gcoeff(ct,j,i) = stoi(Fl_center(cp[t], p, pov2));
    3340             :       }
    3341             :     }
    3342             :   }
    3343             :   else
    3344             :   {
    3345         385 :     const long var = 1;
    3346         385 :     ulong ze = Fl_powu(pgener_Fl(p),(p-1)/expo, p); /* seen as zeta_e^(-1) */
    3347         385 :     GEN vze = Fl_powers(ze, expo-1, p); /* vze[i] = ze^(i-1) */
    3348         385 :     GEN vzeZX = const_vec(p, gen_0);
    3349         385 :     GEN T = polcyclo(expo, var), vT = const_vec(expo,NULL), pro = NULL;
    3350         385 :     long phie = degpol(T), id1 = gel(vjg,1)[1]; /* index of 1_G, always 1 ? */
    3351         385 :     gel(vT, expo) = T;
    3352         385 :     if (f != expo)
    3353             :     {
    3354         147 :       long phif = eulerphiu(f);
    3355         147 :       GEN zf = ZX_rem(pol_xn(expo/f,var), T), zfj = zf;
    3356         147 :       GEN M = cgetg(phif+1, t_MAT);
    3357         147 :       gel(M,1) = col_ei(phie,1);
    3358         518 :       for (j = 2; j <= phif; j++)
    3359             :       {
    3360         371 :         gel(M,j) = RgX_to_RgC(zfj, phie);
    3361         371 :         if (j < phif) zfj = ZX_rem(ZX_mul(zfj, zf), T);
    3362             :       }
    3363         147 :       pro = Qevproj_init(M);
    3364             :     }
    3365         385 :     gel(vzeZX,1) = pol_1(var);
    3366        3416 :     for (i = 2; i <= expo; i++)
    3367             :     {
    3368        3031 :       GEN t = ZX_rem(pol_xn(expo-(i-1), var), T);
    3369        3031 :       if (pro) t = Qevproj_RgX(t, phie, pro);
    3370        3031 :       gel(vzeZX, vze[i]) = t;
    3371             :     }
    3372        3941 :     for(i=1; i<lcl; i++)
    3373             :     { /* loop over characters */
    3374        3556 :       GEN cp = gel(ctp,i), C, cj; /* cp[j] = chi(g_j) mod \P */
    3375        3556 :       long dim = cp[id1];
    3376        3556 :       gel(ct, i) = C = const_col(lcl-1, NULL);
    3377        3556 :       gel(C,operm[1]) = utoi(dim); /* chi(1_G) */
    3378       40978 :       for (j=lcl-1; j > 1; j--)
    3379             :       { /* loop over conjugacy classes, decreasing order: skip 1_G */
    3380       37422 :         long e, jperm = operm[j], o = vord[jperm];
    3381       37422 :         GEN To, jg = gel(vjg,jperm); /* jg[l+1] = class of g^l */
    3382             : 
    3383       37422 :         if (gel(C, jperm)) continue; /* done already */
    3384       35903 :         if (dim == 1) { gel(C, jperm) = gel(vzeZX, cp[jg[2]]); continue; }
    3385         861 :         e = expo / o;
    3386         861 :         cj = chiFT(cp, jg, vze, e, o, p, pov2);
    3387         861 :         To = gel(vT, o); if (!To) To = gel(vT,o) = polcyclo(o, var);
    3388         861 :         gel(C, jperm) = chival(cj, T, To, e, pro, phie);
    3389        3920 :         for (k = 2; k < o; k++)
    3390             :         {
    3391        3059 :           GEN ck = RgX_inflate(cj, k); /* chi(g^k) */
    3392        3059 :           gel(C, jg[k+1]) = chival(ck, T, To, e, pro, phie);
    3393             :         }
    3394             :       }
    3395             :     }
    3396             :   }
    3397         532 :   ct = gen_sort(ct,(void*)cmp_universal,cmp_nodata);
    3398        1736 :   i = 1; while (!vec_isconst(gel(ct,i))) i++;
    3399         532 :   if (i > 1) swap(gel(ct,1), gel(ct,i));
    3400         532 :   return mkvec2(ct, utoipos(f));
    3401             : }
    3402             : GEN
    3403         546 : galoischartable(GEN gal)
    3404             : {
    3405         546 :   pari_sp av = avma;
    3406         546 :   GEN cc = group_to_cc(gal);
    3407         546 :   return gerepilecopy(av, cc_chartable(cc));
    3408             : }
    3409             : 
    3410             : static void
    3411        1491 : checkgaloischar(GEN ch, GEN repr)
    3412             : {
    3413        1491 :   if (gvar(ch) == 0) pari_err_PRIORITY("galoischarpoly",ch,"=",0);
    3414        1491 :   if (!is_vec_t(typ(ch))) pari_err_TYPE("galoischarpoly", ch);
    3415        1491 :   if (lg(repr) != lg(ch)) pari_err_DIM("galoischarpoly");
    3416        1491 : }
    3417             : 
    3418             : static long
    3419        1547 : galoischar_dim(GEN ch)
    3420             : {
    3421        1547 :   pari_sp av = avma;
    3422        1547 :   long d = gtos(simplify_shallow(lift_shallow(gel(ch,1))));
    3423        1547 :   return gc_long(av,d);
    3424             : }
    3425             : 
    3426             : static GEN
    3427       12355 : galoischar_aut_charpoly(GEN cc, GEN ch, GEN p, long d)
    3428             : {
    3429       12355 :   GEN q = p, V = cgetg(d+2, t_POL);
    3430             :   long i;
    3431       12355 :   V[1] = evalsigne(1)|evalvarn(0);
    3432       25970 :   for (i = 1; i <= d; i++)
    3433             :   {
    3434       13615 :     gel(V,i+1) = gel(ch, cc_id(cc,q));
    3435       13615 :     if (i < d) q = perm_mul(q, p);
    3436             :   }
    3437       12355 :   return liftpol_shallow(RgXn_expint(RgX_neg(V),d+1));
    3438             : }
    3439             : 
    3440             : static GEN
    3441        1491 : galoischar_charpoly(GEN cc, GEN ch, long o)
    3442             : {
    3443        1491 :   GEN chm, V, elts = gel(cc,1), repr = gel(cc,3);
    3444        1491 :   long i, d, l = lg(ch), v = gvar(ch);
    3445        1491 :   checkgaloischar(ch, repr);
    3446        1491 :   chm = v < 0 ? ch: gmodulo(ch, polcyclo(o, v));
    3447        1491 :   V = cgetg(l, t_COL); d = galoischar_dim(ch);
    3448       13846 :   for (i = 1; i < l; i++)
    3449       12355 :     gel(V,i) = galoischar_aut_charpoly(cc, chm, gel(elts,repr[i]), d);
    3450        1491 :   return V;
    3451             : }
    3452             : 
    3453             : GEN
    3454        1435 : galoischarpoly(GEN gal, GEN ch, long o)
    3455             : {
    3456        1435 :   pari_sp av = avma;
    3457        1435 :   GEN cc = group_to_cc(gal);
    3458        1435 :   return gerepilecopy(av, galoischar_charpoly(cc, ch, o));
    3459             : }
    3460             : 
    3461             : static GEN
    3462          56 : cc_char_det(GEN cc, GEN ch, long o)
    3463             : {
    3464          56 :   long i, l = lg(ch), d = galoischar_dim(ch);
    3465          56 :   GEN V = galoischar_charpoly(cc, ch, o);
    3466         280 :   for (i = 1; i < l; i++) gel(V,i) = leading_coeff(gel(V,i));
    3467          56 :   return odd(d)? gneg(V): V;
    3468             : }
    3469             : 
    3470             : GEN
    3471          56 : galoischardet(GEN gal, GEN ch, long o)
    3472             : {
    3473          56 :   pari_sp av = avma;
    3474          56 :   GEN cc = group_to_cc(gal);
    3475          56 :   return gerepilecopy(av, cc_char_det(cc, ch, o));
    3476             : }

Generated by: LCOV version 1.13