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