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

Generated by: LCOV version 1.13