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

Generated by: LCOV version 1.16