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

Generated by: LCOV version 1.11