Annotation of OpenXM_contrib/pari-2.2/src/basemath/galconj.c, Revision 1.3
1.3 ! noro 1: /* $Id: galconj.c,v 1.106 2002/09/08 10:47:36 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
1.3 ! noro 16: #include "pari.h"
! 17: extern GEN bezout_lift_fact(GEN T, GEN Tmod, GEN p, long e);
! 18: extern GEN respm(GEN x,GEN y,GEN pm);
! 19: extern GEN ZX_disc_all(GEN,long);
! 20: extern GEN polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom);
! 21: extern GEN znstar_hnf_elts(GEN Z, GEN H);
! 22: extern long ZX_get_prec(GEN x);
! 23:
1.1 noro 24: /*************************************************************************/
25: /** **/
26: /** GALOIS CONJUGATES **/
27: /** **/
28: /*************************************************************************/
1.3 ! noro 29:
1.1 noro 30: GEN
31: galoisconj(GEN nf)
32: {
33: GEN x, y, z;
1.3 ! noro 34: long i, lz, v;
! 35: gpmem_t av = avma;
1.1 noro 36: nf = checknf(nf);
37: x = (GEN) nf[1];
38: v = varn(x);
39: if (v == 0)
40: nf = gsubst(nf, 0, polx[MAXVARN]);
41: else
42: {
43: x = dummycopy(x);
44: setvarn(x, 0);
45: }
46: z = nfroots(nf, x);
47: lz = lg(z);
48: y = cgetg(lz, t_COL);
49: for (i = 1; i < lz; i++)
50: {
51: GEN p1 = lift((GEN) z[i]);
52: setvarn(p1, v);
53: y[i] = (long) p1;
54: }
55: return gerepileupto(av, y);
56: }
57:
58: /* nbmax: maximum number of possible conjugates */
59: GEN
60: galoisconj2pol(GEN x, long nbmax, long prec)
61: {
1.3 ! noro 62: long i, n, v, nbauto;
! 63: gpmem_t av = avma;
1.1 noro 64: GEN y, w, polr, p1, p2;
65: n = degpol(x);
66: if (n <= 0)
67: return cgetg(1, t_VEC);
68: if (gisirreducible(x) == gzero)
69: err(redpoler, "galoisconj2pol");
70: polr = roots(x, prec);
71: p1 = (GEN) polr[1];
72: nbauto = 1;
73: prec = (long) (bit_accuracy(prec) * L2SL10 * 0.75);
74: w = cgetg(n + 2, t_VEC);
75: w[1] = un;
76: for (i = 2; i <= n; i++)
77: w[i] = lmul(p1, (GEN) w[i - 1]);
78: v = varn(x);
79: y = cgetg(nbmax + 1, t_COL);
80: y[1] = (long) polx[v];
81: for (i = 2; i <= n && nbauto < nbmax; i++)
82: {
83: w[n + 1] = polr[i];
84: p1 = lindep2(w, prec);
85: if (signe(p1[n + 1]))
86: {
87: setlg(p1, n + 1);
88: p2 = gdiv(gtopolyrev(p1, v), negi((GEN) p1[n + 1]));
89: if (gdivise(poleval(x, p2), x))
90: {
91: y[++nbauto] = (long) p2;
92: if (DEBUGLEVEL > 1)
93: fprintferr("conjugate %ld: %Z\n", i, y[nbauto]);
94: }
95: }
96: }
97: setlg(y, 1 + nbauto);
98: return gerepileupto(av, gen_sort(y, 0, cmp_pol));
99: }
100:
101: GEN
102: galoisconj2(GEN nf, long nbmax, long prec)
103: {
1.3 ! noro 104: long i, j, n, r1, ru, nbauto;
! 105: gpmem_t av = avma;
1.1 noro 106: GEN x, y, w, polr, p1, p2;
107: if (typ(nf) == t_POL)
108: return galoisconj2pol(nf, nbmax, prec);
109: nf = checknf(nf);
110: x = (GEN) nf[1];
111: n = degpol(x);
112: if (n <= 0)
113: return cgetg(1, t_VEC);
114: r1 = nf_get_r1(nf);
115: p1 = (GEN) nf[6];
116: prec = precision((GEN) p1[1]);
117: /* accuracy in decimal digits */
118: prec = (long) (bit_accuracy(prec) * L2SL10 * 0.75);
119: ru = (n + r1) >> 1;
120: nbauto = 1;
121: polr = cgetg(n + 1, t_VEC);
122: for (i = 1; i <= r1; i++)
123: polr[i] = p1[i];
124: for (j = i; i <= ru; i++)
125: {
126: polr[j++] = p1[i];
127: polr[j++] = lconj((GEN) p1[i]);
128: }
129: p2 = gmael(nf, 5, 1);
130: w = cgetg(n + 2, t_VEC);
131: for (i = 1; i <= n; i++)
132: w[i] = coeff(p2, 1, i);
133: y = cgetg(nbmax + 1, t_COL);
134: y[1] = (long) polx[varn(x)];
135: for (i = 2; i <= n && nbauto < nbmax; i++)
136: {
137: w[n + 1] = polr[i];
138: p1 = lindep2(w, prec);
139: if (signe(p1[n + 1]))
140: {
141: setlg(p1, n + 1);
142: settyp(p1, t_COL);
143: p2 = gdiv(gmul((GEN) nf[7], p1), negi((GEN) p1[n + 1]));
144: if (gdivise(poleval(x, p2), x))
145: {
146: y[++nbauto] = (long) p2;
147: if (DEBUGLEVEL > 1)
148: fprintferr("conjugate %ld: %Z\n", i, y[nbauto]);
149: }
150: }
151: }
152: setlg(y, 1 + nbauto);
153: return gerepileupto(av, gen_sort(y, 0, cmp_pol));
154: }
155: /*************************************************************************/
156: /** **/
157: /** GALOISCONJ4 **/
158: /** **/
159: /** **/
160: /*************************************************************************/
161: /*DEBUGLEVEL:
162: 1: timing
163: 2: outline
164: 4: complete outline
165: 6: detail
166: 7: memory
167: 9: complete detail
168: */
1.3 ! noro 169:
! 170: /* DP = multiple of disc(P) or NULL
! 171: * Return a multiple of the denominator of an algebraic integer (in Q[X]/(P))
! 172: * when expressed in terms of the power basis */
1.1 noro 173: GEN
1.3 ! noro 174: indexpartial(GEN P, GEN DP)
1.1 noro 175: {
1.3 ! noro 176: gpmem_t av = avma;
! 177: long i, nb;
! 178: GEN fa, p1, res = gun, dP;
! 179: pari_timer T;
1.1 noro 180:
1.3 ! noro 181: dP = derivpol(P);
! 182: if(DEBUGLEVEL>=5) (void)TIMER(&T);
! 183: if(!DP) DP = ZX_disc(P);
! 184: DP = mpabs(DP);
! 185: if(DEBUGLEVEL>=5) msgTIMER(&T,"IndexPartial: discriminant");
! 186: fa = auxdecomp(DP, 0);
! 187: if(DEBUGLEVEL>=5) msgTIMER(&T,"IndexPartial: factorization");
! 188: nb = lg(fa[1]);
! 189: for (i = 1; i < nb; i++)
! 190: {
! 191: GEN p=gmael(fa,1,i);
! 192: GEN e=gmael(fa,2,i);
! 193: p1 = powgi(p,shifti(e,-1));
! 194: if ( i==nb-1 )
! 195: {
! 196: if ( mod2(e) && !isprime(p) )
! 197: p1 = mulii(p1,p);
! 198: }
! 199: else if ( cmpis(e,4)>=0 )
! 200: {
! 201: if(DEBUGLEVEL>=5) fprintferr("IndexPartial: factor %Z ",p1);
! 202: p1 = mppgcd(p1, respm(P,dP,p1));
! 203: if(DEBUGLEVEL>=5)
! 204: {
! 205: fprintferr("--> %Z : ",p1);
! 206: msgTIMER(&T,"");
! 207: }
! 208: }
! 209: res=mulii(res,p1);
! 210: }
! 211: return gerepileupto(av,res);
1.1 noro 212: }
213:
1.3 ! noro 214: GEN
! 215: vandermondeinverseprep(GEN L)
1.1 noro 216: {
217: int i, j, n = lg(L);
218: GEN V;
219: V = cgetg(n, t_VEC);
220: for (i = 1; i < n; i++)
221: {
1.3 ! noro 222: gpmem_t ltop=avma;
1.1 noro 223: GEN W=cgetg(n,t_VEC);
224: for (j = 1; j < n; j++)
225: if (i==j)
226: W[j]=un;
227: else
228: W[j]=lsub((GEN)L[i],(GEN)L[j]);
229: V[i]=lpileupto(ltop,divide_conquer_prod(W,&gmul));
230: }
231: return V;
232: }
233:
234: /* Calcule l'inverse de la matrice de van der Monde de T multiplie par den */
235: GEN
236: vandermondeinverse(GEN L, GEN T, GEN den, GEN prep)
237: {
1.3 ! noro 238: gpmem_t ltop = avma;
! 239: int i, n = lg(L)-1;
1.1 noro 240: long x = varn(T);
241: GEN M, P;
242: if (!prep)
1.3 ! noro 243: prep = vandermondeinverseprep(L);
! 244: M = cgetg(n+1, t_MAT);
! 245: for (i = 1; i <= n; i++)
1.1 noro 246: {
247: P = gdiv(gdeuc(T, gsub(polx[x], (GEN) L[i])), (GEN) prep[i]);
1.3 ! noro 248: M[i] = (long)pol_to_vec(P,n);
1.1 noro 249: }
1.3 ! noro 250: return gerepileupto(ltop, gmul(den, M));
1.1 noro 251: }
252:
253: /* Calcule les bornes sur les coefficients a chercher */
254: struct galois_borne
255: {
256: GEN l;
257: long valsol;
258: long valabs;
259: GEN bornesol;
260: GEN ladicsol;
261: GEN ladicabs;
262: GEN lbornesol;
263: };
264:
265:
266: GEN
1.3 ! noro 267: initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis)
1.1 noro 268: {
1.3 ! noro 269: int i, n = degpol(T);
! 270: GEN L, z, prep, den;
! 271: pari_timer ti;
! 272:
! 273: if (DEBUGLEVEL>=4) (void)TIMER(&ti);
1.1 noro 274: L = roots(T, prec);
1.3 ! noro 275: if (DEBUGLEVEL>=4) msgTIMER(&ti,"roots");
1.1 noro 276: for (i = 1; i <= n; i++)
277: {
278: z = (GEN) L[i];
1.3 ! noro 279: if (signe(z[2])) break;
1.1 noro 280: L[i] = z[1];
281: }
1.3 ! noro 282: if (DEBUGLEVEL>=4) (void)TIMER(&ti);
! 283: prep = vandermondeinverseprep(L);
1.1 noro 284: if (!dn)
285: {
1.3 ! noro 286: GEN dis, res = divide_conquer_prod(gabs(prep,prec), mpmul);
1.1 noro 287: disable_dbg(0);
1.3 ! noro 288: dis = ZX_disc_all(T, 1+logint(res,gdeux,NULL));
1.1 noro 289: disable_dbg(-1);
1.3 ! noro 290: den = indexpartial(T,dis);
! 291: if (ptdis) *ptdis = dis;
1.1 noro 292: }
293: else
1.3 ! noro 294: den = dn;
! 295: if (ptprep) *ptprep = prep;
! 296: *ptL = L; return den;
! 297: }
! 298:
! 299: /* ||| M ||| with respect to || x ||_oo. Assume M square t_MAT */
! 300: GEN
! 301: matrixnorm(GEN M, long prec)
! 302: {
! 303: long i,j, n = lg(M);
! 304: GEN B = realzero(prec);
! 305:
! 306: for (i = 1; i < n; i++)
1.1 noro 307: {
1.3 ! noro 308: GEN z = gabs(gcoeff(M,i,1), prec);
! 309: for (j = 2; j < n; j++)
1.1 noro 310: z = gadd(z, gabs(gcoeff(M,i,j), prec));
1.3 ! noro 311: if (gcmp(z, B) > 0) B = z;
1.1 noro 312: }
1.3 ! noro 313: return B;
! 314: }
! 315:
! 316: /* L a t_VEC/t_COL, return ||L||_oo */
! 317: GEN
! 318: supnorm(GEN L, long prec)
! 319: {
! 320: long i, n = lg(L);
! 321: GEN z, B;
! 322:
! 323: if (n == 1) return realzero(prec);
! 324: B = gabs((GEN)L[1], prec);
! 325: for (i = 2; i < n; i++)
1.1 noro 326: {
1.3 ! noro 327: z = gabs((GEN)L[i], prec);
! 328: if (gcmp(z, B) > 0) B = z;
1.1 noro 329: }
1.3 ! noro 330: return B;
! 331: }
! 332:
! 333: GEN
! 334: galoisborne(GEN T, GEN dn, struct galois_borne *gb, long ppp)
! 335: {
! 336: gpmem_t ltop = avma, av2;
! 337: GEN borne, borneroots, borneabs;
! 338: int n;
! 339: GEN L, M, prep, den;
! 340: long prec;
! 341: pari_timer ti;
! 342:
! 343: prec = ZX_get_prec(T);
! 344: den = initgaloisborne(T,dn,prec, &L,&prep,NULL);
! 345: if (!dn) den = gclone(den);
! 346: if (DEBUGLEVEL>=4) TIMERstart(&ti);
! 347: M = vandermondeinverse(L, gmul(T, realun(prec)), den, prep);
! 348: if (DEBUGLEVEL>=4) msgTIMER(&ti,"vandermondeinverse");
! 349: borne = matrixnorm(M, prec);
! 350: borneroots = supnorm(L, prec);
! 351: n = degpol(T);
1.1 noro 352: borneabs = addsr(1, gmulsg(n, gpowgs(borneroots, n/ppp)));
353: /*if (ppp == 1)
354: borneabs = addsr(1, gmulsg(n, gpowgs(borneabs, 2)));*/
355: borneroots = addsr(1, gmul(borne, borneroots));
356: av2 = avma;
357: /*We use d-1 test, so we must overlift to 2^BITS_IN_LONG*/
1.3 ! noro 358: gb->valsol = logint(gmul2n(borneroots,2+BITS_IN_LONG), gb->l,NULL);
! 359: gb->valabs = logint(gmul2n(borneabs,2), gb->l,NULL);
! 360: gb->valabs = max(gb->valsol, gb->valabs);
1.1 noro 361: if (DEBUGLEVEL >= 4)
362: fprintferr("GaloisConj:val1=%ld val2=%ld\n", gb->valsol, gb->valabs);
363: avma = av2;
1.3 ! noro 364: gb->bornesol = gerepileupto(ltop, ceil_safe(mulrs(borneroots,2)));
1.1 noro 365: if (DEBUGLEVEL >= 9)
366: fprintferr("GaloisConj: Bound %Z\n",borneroots);
367: gb->ladicsol = gpowgs(gb->l, gb->valsol);
368: gb->ladicabs = gpowgs(gb->l, gb->valabs);
369: gb->lbornesol = subii(gb->ladicsol,gb->bornesol);
1.3 ! noro 370: if (!dn) { dn = icopy(den); gunclone(den); }
1.1 noro 371: return dn;
372: }
373:
374: struct galois_lift
375: {
376: GEN T;
377: GEN den;
378: GEN p;
379: GEN L;
380: GEN Lden;
381: long e;
382: GEN Q;
383: GEN TQ;
384: struct galois_borne *gb;
385: };
386: /* Initialise la structure galois_lift */
387: GEN makeLden(GEN L,GEN den, struct galois_borne *gb)
388: {
1.3 ! noro 389: gpmem_t ltop=avma;
1.1 noro 390: long i,l=lg(L);
391: GEN Lden;
392: Lden=cgetg(l,t_VEC);
393: for (i=1;i<l;i++)
394: Lden[i]=lmulii((GEN)L[i],den);
395: for (i=1;i<l;i++)
396: Lden[i]=lmodii((GEN)Lden[i],gb->ladicsol);
397: return gerepileupto(ltop,Lden);
398: }
399: void
400: initlift(GEN T, GEN den, GEN p, GEN L, GEN Lden, struct galois_borne *gb, struct galois_lift *gl)
401: {
1.3 ! noro 402: gpmem_t ltop, lbot;
1.1 noro 403: gl->gb=gb;
404: gl->T = T;
405: gl->den = den;
406: gl->p = p;
407: gl->L = L;
408: gl->Lden = Lden;
409: ltop = avma;
1.3 ! noro 410: gl->e = logint(gmul2n(gb->bornesol, 2+BITS_IN_LONG),p,NULL);
1.1 noro 411: gl->e = max(2,gl->e);
412: lbot = avma;
413: gl->Q = gpowgs(p, gl->e);
414: gl->Q = gerepile(ltop, lbot, gl->Q);
415: gl->TQ = FpX_red(T,gl->Q);
416: }
417:
418: /*
419: * Verifie que f est une solution presque surement et calcule sa permutation
420: */
421: static int
422: poltopermtest(GEN f, struct galois_lift *gl, GEN pf)
423: {
1.3 ! noro 424: gpmem_t ltop;
1.1 noro 425: GEN fx, fp;
426: long i, j,ll;
427: for (i = 2; i< lgef(f); i++)
428: if (cmpii((GEN)f[i],gl->gb->bornesol)>0
429: && cmpii((GEN)f[i],gl->gb->lbornesol)<0)
430: {
431: if (DEBUGLEVEL>=4)
432: fprintferr("GaloisConj: Solution too large, discard it.\n");
1.3 ! noro 433: if (DEBUGLEVEL>=8)
! 434: fprintferr("f=%Z\n borne=%Z\n l-borne=%Z\n",f,gl->gb->bornesol,gl->gb->lbornesol);
1.1 noro 435: return 0;
436: }
437: ll=lg(gl->L);
438: fp = cgetg(ll, t_VECSMALL);
439: ltop = avma;
440: for (i = 1; i < ll; i++)
441: fp[i] = 1;
442: for (i = 1; i < ll; i++)
443: {
444: fx = FpX_eval(f, (GEN) gl->L[i],gl->gb->ladicsol);
445: for (j = 1; j < ll; j++)
446: {
447: if (fp[j] && egalii(fx, (GEN) gl->Lden[j]))
448: {
449: pf[i] = j;
450: fp[j] = 0;
451: break;
452: }
453: }
454: if (j == ll)
455: return 0;
456: avma = ltop;
457: }
458: return 1;
459: }
460:
461: /*
462: * Soit P un polynome de \ZZ[X] , p un nombre premier , S\in\FF_p[X]/(Q) tel
463: * que P(S)=0 [p,Q] Relever S en S_0 tel que P(S_0)=0 [p^e,Q]
464: * Unclean stack.
465: */
466: static long monoratlift(GEN S, GEN q, GEN qm1old,struct galois_lift *gl, GEN frob)
467: {
1.3 ! noro 468: gpmem_t ltop=avma;
1.1 noro 469: GEN tlift = polratlift(S,q,qm1old,qm1old,gl->den);
470: if (tlift)
471: {
472: if(DEBUGLEVEL>=4)
473: fprintferr("MonomorphismLift: trying early solution %Z\n",tlift);
474: /*Rationals coefficients*/
475: tlift=lift(gmul(tlift,gmodulcp(gl->den,gl->gb->ladicsol)));
476: if (poltopermtest(tlift, gl, frob))
477: {
478: if(DEBUGLEVEL>=4)
479: fprintferr("MonomorphismLift: true early solution.\n");
480: avma=ltop;
481: return 1;
482: }
483: if(DEBUGLEVEL>=4)
484: fprintferr("MonomorphismLift: false early solution.\n");
485: }
486: avma=ltop;
487: return 0;
488: }
489: GEN
490: monomorphismratlift(GEN P, GEN S, struct galois_lift *gl, GEN frob)
491: {
1.3 ! noro 492: gpmem_t ltop, lbot;
1.1 noro 493: long rt;
494: GEN Q=gl->T, p=gl->p;
495: long e=gl->e, level=1;
496: long x;
497: GEN q, qold, qm1, qm1old;
498: GEN W, Pr, Qr, Sr, Wr = gzero, Prold, Qrold, Spow;
499: long i,nb,mask;
500: GEN *gptr[2];
1.3 ! noro 501: if (DEBUGLEVEL == 1) (void)timer2();
1.1 noro 502: x = varn(P);
503: rt = brent_kung_optpow(degpol(Q),1);
504: q = p; qm1 = gun; /*during the run, we have p*qm1=q*/
505: nb=hensel_lift_accel(e, &mask);
506: Pr = FpX_red(P,q);
507: Qr = (P==Q)?Pr:FpX_red(Q, q);/*A little speed up for automorphismlift*/
508: W=FpX_FpXQ_compo(deriv(Pr, x),S,Qr,q);
509: W=FpXQ_inv(W,Qr,q);
510: qold = p; qm1old=gun;
511: Prold = Pr;
512: Qrold = Qr;
513: gptr[0] = &S;
514: gptr[1] = &Wr;
515: for (i=0; i<nb;i++)
516: {
517: if (DEBUGLEVEL>=2)
518: {
519: level=(level<<1)-((mask>>i)&1);
1.3 ! noro 520: (void)timer2();
1.1 noro 521: }
522: qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
523: q = mulii(qm1, p);
524: Pr = FpX_red(P, q);
525: Qr = (P==Q)?Pr:FpX_red(Q, q);/*A little speed up for automorphismlift*/
526: ltop = avma;
527: Sr = S;
528: Spow = FpXQ_powers(Sr, rt, Qr, q);
529:
530: if (i)
531: {
532: W = FpXQ_mul(Wr, FpX_FpXQV_compo(deriv(Pr,-1),FpXV_red(Spow,qold),Qrold,qold), Qrold, qold);
533: W = FpX_neg(W, qold);
534: W = FpX_Fp_add(W, gdeux, qold);
535: W = FpXQ_mul(Wr, W, Qrold, qold);
536: }
537: Wr = W;
538: S = FpXQ_mul(Wr, FpX_FpXQV_compo(Pr, Spow, Qr, q),Qr,q);
539: S = FpX_sub(Sr, S, NULL);
540: lbot = avma;
541: Wr = gcopy(Wr);
542: S = FpX_red(S, q);
543: gerepilemanysp(ltop, lbot, gptr, 2);
544: if (i && i<nb-1 && frob && monoratlift(S,q,qm1old,gl,frob))
545: return NULL;
546: qold = q; qm1old=qm1; Prold = Pr; Qrold = Qr;
547: if (DEBUGLEVEL >= 2)
548: msgtimer("MonomorphismLift: lift to prec %d",level);
549: }
550: if (DEBUGLEVEL == 1)
551: msgtimer("monomorphismlift()");
552: return S;
553: }
554: /*
555: * Soit T un polynome de \ZZ[X] , p un nombre premier , S\in\FF_p[X]/(T) tel
556: * que T(S)=0 [p,T] Relever S en S_0 tel que T(S_0)=0 [T,p^e]
557: * Unclean stack.
558: */
559: GEN
560: automorphismlift(GEN S, struct galois_lift *gl, GEN frob)
561: {
562: return monomorphismratlift(gl->T, S, gl, frob);
563: }
564: GEN
565: monomorphismlift(GEN P, GEN S, GEN Q, GEN p, long e)
566: {
567: struct galois_lift gl;
568: gl.T=Q;
569: gl.p=p;
570: gl.e=e;
571: return monomorphismratlift(P,S,&gl,NULL);
572: }
573:
574: struct galois_testlift
575: {
576: long n;
577: long f;
578: long g;
579: GEN bezoutcoeff;
580: GEN pauto;
581: GEN C;
582: GEN Cd;
583: };
584: static GEN
585: galoisdolift(struct galois_lift *gl, GEN frob)
586: {
1.3 ! noro 587: gpmem_t ltop=avma;
1.1 noro 588: long v = varn(gl->T);
589: GEN Tp = FpX_red(gl->T, gl->p);
590: GEN S = FpXQ_pow(polx[v],gl->p, Tp,gl->p);
591: GEN plift = automorphismlift(S, gl, frob);
592: return gerepileupto(ltop,plift);
593: }
594:
595: static void
596: inittestlift( GEN plift, GEN Tmod, struct galois_lift *gl,
597: struct galois_testlift *gt)
598: {
599: long v = varn(gl->T);
600: gt->n = lg(gl->L) - 1;
601: gt->g = lg(Tmod) - 1;
602: gt->f = gt->n / gt->g;
603: gt->bezoutcoeff = bezout_lift_fact(gl->T, Tmod, gl->p, gl->e);
604: gt->pauto = cgetg(gt->f + 1, t_VEC);
605: gt->pauto[1] = (long) polx[v];
606: gt->pauto[2] = (long) gcopy(plift);
607: if (gt->f > 2)
608: {
1.3 ! noro 609: gpmem_t ltop=avma;
1.1 noro 610: int i;
611: long nautpow=brent_kung_optpow(gt->n-1,gt->f-2);
612: GEN autpow;
1.3 ! noro 613: if (DEBUGLEVEL >= 1) (void)timer2();
1.1 noro 614: autpow = FpXQ_powers(plift,nautpow,gl->TQ,gl->Q);
615: for (i = 3; i <= gt->f; i++)
616: gt->pauto[i] = (long) FpX_FpXQV_compo((GEN)gt->pauto[i-1],autpow,gl->TQ,gl->Q);
617: /*Somewhat paranoid with memory, but this function use a lot of stack.*/
618: gt->pauto=gerepileupto(ltop, gt->pauto);
619: if (DEBUGLEVEL >= 1) msgtimer("frobenius power");
620: }
621: }
622:
623: long intheadlong(GEN x, GEN mod)
624: {
1.3 ! noro 625: gpmem_t ltop=avma;
1.1 noro 626: GEN r;
627: int s;
628: long res;
629: r=divii(shifti(x,BITS_IN_LONG),mod);
630: s=signe(r);
631: res=s?s>0?r[2]:-r[2]:0;
632: avma=ltop;
633: return res;
634: }
635:
636: GEN matheadlong(GEN W, GEN mod)
637: {
638: long i,j;
639: GEN V=cgetg(lg(W),t_VEC);
640: for(i=1;i<lg(W);i++)
641: {
642: V[i]=lgetg(lg(W[i]),t_VECSMALL);
643: for(j=1;j<lg(W[i]);j++)
644: mael(V,i,j)=intheadlong(gmael(W,i,j),mod);
645: }
646: return V;
647: }
648:
649: long polheadlong(GEN P, long n, GEN mod)
650: {
651: return (lgef(P)>n+2)?intheadlong((GEN)P[n+2],mod):0;
652: }
653: /*
654: *
655: */
656: long
657: frobeniusliftall(GEN sg, long el, GEN *psi, struct galois_lift *gl,
658: struct galois_testlift *gt, GEN frob)
659: {
1.3 ! noro 660: gpmem_t av, ltop2, ltop = avma;
1.1 noro 661: long d, z, m, c, n, ord;
662: int i, j, k;
663: GEN pf, u, v;
1.3 ! noro 664: GEN C, Cd, sgi, cache;
! 665: long Z, c_idx=gt->g-1;
! 666: long stop=0,hop=0;
! 667: GEN NN,NQ,NR;
! 668: long N1,N2,R1,Ni;
1.1 noro 669: m = gt->g;
670: ord = gt->f;
671: n = lg(gl->L) - 1;
672: c = lg(sg) - 1;
673: d = m / c;
674: pf = cgetg(m, t_VECSMALL);
675: *psi = pf;
676: ltop2 = avma;
677: NN = gdiv(mpfact(m), gmul(stoi(c), gpowgs(mpfact(d), c)));
678: if (DEBUGLEVEL >= 4)
679: fprintferr("GaloisConj:I will try %Z permutations\n", NN);
680: N1=10000000;
681: NQ=dvmdis(NN,N1,&NR);
682: if (cmpis(NQ,1000000000)>0)
683: {
1.3 ! noro 684: err(warner,"Combinatorics too hard : would need %Z tests!\n"
! 685: "I will skip it,but it may induce galoisinit to loop",NN);
1.1 noro 686: avma = ltop;
687: *psi = NULL;
688: return 0;
689: }
690: N2=itos(NQ); R1=itos(NR); if(!N2) N1=R1;
691: if (DEBUGLEVEL>=4)
692: {
693: stop=N1/20;
1.3 ! noro 694: (void)timer2();
1.1 noro 695: }
696: avma = ltop2;
697: C=gt->C;
698: Cd=gt->Cd;
699: v = FpX_Fp_mul(FpXQ_mul((GEN)gt->pauto[1+el%ord], (GEN)
700: gt->bezoutcoeff[m],gl->TQ,gl->Q),gl->den,gl->Q);
1.3 ! noro 701: sgi=cgetg(lg(sg),t_VECSMALL);
! 702: for(i=1;i<lg(sgi);i++)
! 703: sgi[i]=(el*sg[i])%ord + 1;
! 704: cache=cgetg(m+1,t_VECSMALL);
! 705: cache[m]=polheadlong(v,1,gl->Q);
! 706: Z=polheadlong(v,2,gl->Q);
1.1 noro 707: for (i = 1; i < m; i++)
708: pf[i] = 1 + i / d;
709: av = avma;
710: for (Ni = 0, i = 0; ;i++)
711: {
1.3 ! noro 712: for (j = c_idx ; j > 0; j--)
1.1 noro 713: {
714: long h;
1.3 ! noro 715: h=sgi[pf[j]];
1.1 noro 716: if (!mael(C,h,j))
717: {
1.3 ! noro 718: gpmem_t av3=avma;
1.1 noro 719: GEN r;
720: r=FpX_Fp_mul(FpXQ_mul((GEN) gt->pauto[h],
721: (GEN) gt->bezoutcoeff[j],gl->TQ,gl->Q),gl->den,gl->Q);
722: mael(C,h,j) = lclone(r);
1.3 ! noro 723: mael(Cd,h,j) = polheadlong(r,1,gl->Q);
1.1 noro 724: avma=av3;
725: }
1.3 ! noro 726: cache[j]=cache[j+1]+mael(Cd,h,j);
1.1 noro 727: }
1.3 ! noro 728: if (labs(cache[1])<=n )
1.1 noro 729: {
1.3 ! noro 730: long ZZ=Z;
1.1 noro 731: for (j = 1; j < m; j++)
1.3 ! noro 732: ZZ += polheadlong(gmael(C,sgi[pf[j]],j),2,gl->Q);
! 733: if (labs(ZZ)<=n )
1.1 noro 734: {
735: u = v;
736: for (j = 1; j < m; j++)
1.3 ! noro 737: u = FpX_add(u, gmael(C,sgi[pf[j]],j),NULL);
1.1 noro 738: u = FpX_center(FpX_red(u, gl->Q), gl->Q);
739: if (poltopermtest(u, gl, frob))
740: {
741: if (DEBUGLEVEL >= 4 )
1.3 ! noro 742: {
1.1 noro 743: msgtimer("");
1.3 ! noro 744: fprintferr("GaloisConj: %d hops on %Z tests\n",hop,addis(mulss(Ni,N1),i));
! 745: }
1.1 noro 746: avma = ltop2;
747: return 1;
748: }
1.3 ! noro 749: else if (DEBUGLEVEL >= 4 )
! 750: fprintferr("M");
1.1 noro 751: }
1.3 ! noro 752: else hop++;
1.1 noro 753: }
754: if (DEBUGLEVEL >= 4 && i==stop)
755: {
756: stop+=N1/20;
757: msgtimer("GaloisConj:Testing %Z", addis(mulss(Ni,N1),i));
758: }
759: avma = av;
760: if (i == N1 - 1)
761: {
762: if (Ni==N2-1)
763: N1=R1;
764: if (Ni==N2)
765: break;
766: Ni++;
767: i=0;
768: if (DEBUGLEVEL>=4)
769: {
770: stop=N1/20;
1.3 ! noro 771: (void)timer2();
1.1 noro 772: }
773: }
774: for (j = 2; j < m && pf[j - 1] >= pf[j]; j++);
775: for (k = 1; k < j - k && pf[k] != pf[j - k]; k++)
776: {
777: z = pf[k];
778: pf[k] = pf[j - k];
779: pf[j - k] = z;
780: }
781: for (k = j - 1; pf[k] >= pf[j]; k--);
782: z = pf[j];
783: pf[j] = pf[k];
784: pf[k] = z;
1.3 ! noro 785: c_idx=j;
1.1 noro 786: }
1.3 ! noro 787: if (DEBUGLEVEL>=4)
! 788: fprintferr("GaloisConj: not found, %d hops \n",hop);
1.1 noro 789: avma = ltop;
790: *psi = NULL;
791: return 0;
792: }
793:
794: /*
795: * alloue une ardoise pour n entiers de longueurs pour les test de
796: * permutation
797: */
798: GEN
799: alloue_ardoise(long n, long s)
800: {
801: GEN ar;
802: long i;
803: ar = cgetg(n + 1, t_VEC);
804: for (i = 1; i <= n; i++)
805: ar[i] = lgetg(s, t_INT);
806: return ar;
807: }
808:
809: /*
810: * structure contenant toutes les données pour le tests des permutations:
811: *
812: * ordre :ordre des tests pour verifie_test ordre[lg(ordre)]: numero du test
813: * principal borne : borne sur les coefficients a trouver ladic: modulo
814: * l-adique des racines lborne:ladic-borne TM:vecteur des ligne de M
815: * PV:vecteur des clones des matrices de test (Vmatrix) (ou NULL si non
816: * calcule) L,M comme d'habitude (voir plus bas)
817: */
818: struct galois_test
819: {
820: GEN ordre;
821: GEN borne, lborne, ladic;
822: GEN PV, TM;
823: GEN L, M;
824: };
825: /* Calcule la matrice de tests correspondant a la n-ieme ligne de V */
826: static GEN
827: Vmatrix(long n, struct galois_test *td)
828: {
1.3 ! noro 829: gpmem_t lbot, ltop = avma;
1.1 noro 830: GEN V;
831: long i;
832: V = cgetg(lg(td->L), t_VEC);
833: for (i = 1; i < lg(V); i++)
834: V[i] = mael(td->M,i,n);
835: V = gmul(td->L, V);
836: lbot = avma;
837: V = gmod(V, td->ladic);
838: return gerepile(ltop, lbot, V);
839: }
840:
841: /*
842: * Initialise la structure galois_test
843: */
844: static void
845: inittest(GEN L, GEN M, GEN borne, GEN ladic, struct galois_test *td)
846: {
1.3 ! noro 847: gpmem_t ltop;
1.1 noro 848: long i;
849: int n = lg(L) - 1;
850: if (DEBUGLEVEL >= 8)
851: fprintferr("GaloisConj:Entree Init Test\n");
852: td->ordre = cgetg(n + 1, t_VECSMALL);
853: for (i = 1; i <= n - 2; i++)
854: td->ordre[i] = i + 2;
855: for (; i <= n; i++)
856: td->ordre[i] = i - n + 2;
857: td->borne = borne;ltop = avma;
858: td->lborne = subii(ladic, borne);
859: td->ladic = ladic;
860: td->L = L;
861: td->M = M;
862: td->PV = cgetg(n + 1, t_VECSMALL); /* peut-etre t_VEC est correct ici */
863: for (i = 1; i <= n; i++)
864: td->PV[i] = 0;
865: ltop = avma;
866: td->PV[td->ordre[n]] = (long) gclone(Vmatrix(td->ordre[n], td));
867: avma = ltop;
1.3 ! noro 868: td->TM = gtrans_i(M);
1.1 noro 869: settyp(td->TM, t_VEC);
870: for (i = 1; i < lg(td->TM); i++)
871: settyp(td->TM[i], t_VEC);
872: if (DEBUGLEVEL >= 8)
873: fprintferr("GaloisConj:Sortie Init Test\n");
874: }
875:
876: /*
877: * liberer les clones de la structure galois_test
878: *
879: * Reservé a l'accreditation ultra-violet:Liberez les clones! Paranoia(tm)
880: */
881: static void
882: freetest(struct galois_test *td)
883: {
884: int i;
885: for (i = 1; i < lg(td->PV); i++)
886: if (td->PV[i])
887: {
888: gunclone((GEN) td->PV[i]);
889: td->PV[i] = 0;
890: }
891: }
892:
893: /*
894: * Test si le nombre padique P est proche d'un entier inferieur a td->borne
895: * en valeur absolue.
896: */
897: static long
898: padicisint(GEN P, struct galois_test *td)
899: {
900: long r;
1.3 ! noro 901: gpmem_t ltop = avma;
1.1 noro 902: GEN U;
903: /*if (typ(P) != t_INT) err(typeer, "padicisint");*/
904: U = modii(P, td->ladic);
905: r = gcmp(U, (GEN) td->borne) <= 0 || gcmp(U, (GEN) td->lborne) >= 0;
906: avma = ltop;
907: return r;
908: }
909:
910: /*
911: * Verifie si pf est une vrai solution et non pas un "hop"
912: */
913: static long
914: verifietest(GEN pf, struct galois_test *td)
915: {
1.3 ! noro 916: gpmem_t av = avma;
1.1 noro 917: GEN P, V;
918: int i, j;
919: int n = lg(td->L) - 1;
920: if (DEBUGLEVEL >= 8)
921: fprintferr("GaloisConj:Entree Verifie Test\n");
1.3 ! noro 922: P = perm_mul(td->L, pf);
1.1 noro 923: for (i = 1; i < n; i++)
924: {
925: long ord;
926: GEN PW;
927: ord = td->ordre[i];
928: PW = (GEN) td->PV[ord];
929: if (PW)
930: {
931: V = ((GEN **) PW)[1][pf[1]];
932: for (j = 2; j <= n; j++)
933: V = addii(V, ((GEN **) PW)[j][pf[j]]);
934: }
935: else
936: V = centermod(gmul((GEN) td->TM[ord], P),td->ladic);
937: if (!padicisint(V, td))
938: break;
939: }
940: if (i == n)
941: {
942: if (DEBUGLEVEL >= 8)
943: fprintferr("GaloisConj:Sortie Verifie Test:1\n");
944: avma = av;
945: return 1;
946: }
947: if (!td->PV[td->ordre[i]])
948: {
949: td->PV[td->ordre[i]] = (long) gclone(Vmatrix(td->ordre[i], td));
950: if (DEBUGLEVEL >= 4)
951: fprintferr("M");
952: }
953: if (DEBUGLEVEL >= 4)
954: fprintferr("%d.", i);
955: if (i > 1)
956: {
957: long z;
958: z = td->ordre[i];
959: for (j = i; j > 1; j--)
960: td->ordre[j] = td->ordre[j - 1];
961: td->ordre[1] = z;
962: if (DEBUGLEVEL >= 8)
963: fprintferr("%Z", td->ordre);
964: }
965: if (DEBUGLEVEL >= 8)
966: fprintferr("GaloisConj:Sortie Verifie Test:0\n");
967: avma = av;
968: return 0;
969: }
970: /*Compute a*b/c when a*b will overflow*/
971: static long muldiv(long a,long b,long c)
972: {
973: return (long)((double)a*(double)b/c);
974: }
975:
976: /*
977: * F est la decomposition en cycle de sigma, B est la decomposition en cycle
978: * de cl(tau) Teste toutes les permutations pf pouvant correspondre a tau tel
979: * que : tau*sigma*tau^-1=sigma^s et tau^d=sigma^t ou d est l'ordre de
980: * cl(tau) --------- x: vecteur des choix y: vecteur des mises a jour G:
981: * vecteur d'acces a F linéaire
982: */
983: GEN
984: testpermutation(GEN F, GEN B, long s, long t, long cut,
985: struct galois_test *td)
986: {
1.3 ! noro 987: gpmem_t av, avm = avma;
1.1 noro 988: int a, b, c, d, n;
989: GEN pf, x, ar, y, *G;
990: int p1, p2, p3, p4, p5, p6;
991: long l1, l2, N1, N2, R1;
992: long i, j, cx, hop = 0, start = 0;
993: GEN W, NN, NQ, NR;
994: long V;
1.3 ! noro 995: if (DEBUGLEVEL >= 1) (void)timer2();
1.1 noro 996: a = lg(F) - 1;
997: b = lg(F[1]) - 1;
998: c = lg(B) - 1;
999: d = lg(B[1]) - 1;
1000: n = a * b;
1001: s = (b + s) % b;
1002: pf = cgetg(n + 1, t_VECSMALL);
1003: av = avma;
1004: ar = cgetg(a + 1, t_VECSMALL);
1005: x = cgetg(a + 1, t_VECSMALL);
1006: y = cgetg(a + 1, t_VECSMALL);
1007: G = (GEN *) cgetg(a + 1, t_VECSMALL); /* Don't worry */
1008: W = matheadlong((GEN) td->PV[td->ordre[n]],td->ladic);
1009: for (cx = 1, i = 1, j = 1; cx <= a; cx++, i++)
1010: {
1011: x[cx] = (i != d) ? 0 : t;
1012: y[cx] = 1;
1013: G[cx] = (GEN) F[((long **) B)[j][i]]; /* Be happy */
1014: if (i == d)
1015: {
1016: i = 0;
1017: j++;
1018: }
1019: } /* Be happy now! */
1020: NN = divis(gpowgs(stoi(b), c * (d - 1)),cut);
1021: if (DEBUGLEVEL >= 4)
1022: fprintferr("GaloisConj:I will try %Z permutations\n", NN);
1023: N1=100000;
1024: NQ=dvmdis(NN,N1,&NR);
1025: if (cmpis(NQ,1000000000)>0)
1026: {
1027: avma=avm;
1028: err(warner,"Combinatorics too hard : would need %Z tests!\n I'll skip it but you will get a partial result...",NN);
1.3 ! noro 1029: return perm_identity(n);
1.1 noro 1030: }
1031: N2=itos(NQ); R1=itos(NR);
1032: for (l2 = 0; l2 <= N2; l2++)
1033: {
1034: long nbiter = (l2<N2) ? N1: R1;
1035: if (DEBUGLEVEL >= 2 && N2)
1036: fprintferr("%d%% ", muldiv(l2,100,N2));
1037: for (l1 = 0; l1 < nbiter; l1++)
1038: {
1039: if (start)
1040: {
1041: for (i = 1, j = d; i < a;)
1042: {
1043: y[i] = 1;
1044: if ((++(x[i])) != b)
1045: break;
1046: x[i++] = 0;
1047: if (i == j)
1048: {
1049: y[i++] = 1;
1050: j += d;
1051: }
1052: }
1053: y[i + 1] = 1;
1054: }
1055: else start=1;
1056: for (p1 = 1, p5 = d; p1 <= a; p1++, p5++)
1057: if (y[p1])
1058: {
1059: if (p5 == d)
1060: {
1061: p5 = 0;
1062: p4 = 0;
1063: }
1064: else
1065: p4 = x[p1 - 1];
1066: if (p5 == d - 1)
1067: p6 = p1 + 1 - d;
1068: else
1069: p6 = p1 + 1;
1070: y[p1] = 0;
1071: V = 0;
1072: for (p2 = 1 + p4, p3 = 1 + x[p1]; p2 <= b; p2++)
1073: {
1074: V += mael(W,G[p6][p3],G[p1][p2]);
1075: p3 += s;
1076: if (p3 > b)
1077: p3 -= b;
1078: }
1079: p3 = 1 + x[p1] - s;
1080: if (p3 <= 0)
1081: p3 += b;
1082: for (p2 = p4; p2 >= 1; p2--)
1083: {
1084: V += mael(W,G[p6][p3],G[p1][p2]);
1085: p3 -= s;
1086: if (p3 <= 0)
1087: p3 += b;
1088: }
1089: ar[p1]=V;
1090: }
1091: V = 0;
1092: for (p1 = 1; p1 <= a; p1++)
1093: V += ar[p1];
1094: if (labs(V)<=n)
1095: {
1096: for (p1 = 1, p5 = d; p1 <= a; p1++, p5++)
1097: {
1098: if (p5 == d)
1099: {
1100: p5 = 0;
1101: p4 = 0;
1102: }
1103: else
1104: p4 = x[p1 - 1];
1105: if (p5 == d - 1)
1106: p6 = p1 + 1 - d;
1107: else
1108: p6 = p1 + 1;
1109: for (p2 = 1 + p4, p3 = 1 + x[p1]; p2 <= b; p2++)
1110: {
1111: pf[G[p1][p2]] = G[p6][p3];
1112: p3 += s;
1113: if (p3 > b)
1114: p3 -= b;
1115: }
1116: p3 = 1 + x[p1] - s;
1117: if (p3 <= 0)
1118: p3 += b;
1119: for (p2 = p4; p2 >= 1; p2--)
1120: {
1121: pf[G[p1][p2]] = G[p6][p3];
1122: p3 -= s;
1123: if (p3 <= 0)
1124: p3 += b;
1125: }
1126: }
1127: if (verifietest(pf, td))
1128: {
1129: if (DEBUGLEVEL >= 1)
1130: {
1131: GEN nb=addis(mulss(l2,N1),l1);
1132: msgtimer("testpermutation(%Z)", nb);
1133: if (DEBUGLEVEL >= 2 && hop)
1134: fprintferr("GaloisConj:%d hop sur %Z iterations\n", hop, nb);
1135: }
1136: avma = av;
1137: return pf;
1138: }
1139: else
1140: hop++;
1141: }
1142: }
1143: }
1144: if (DEBUGLEVEL >= 1)
1145: {
1146: msgtimer("testpermutation(%Z)", NN);
1147: if (DEBUGLEVEL >= 2 && hop)
1148: fprintferr("GaloisConj:%d hop sur %Z iterations\n", hop, NN);
1149: }
1150: avma = avm;
1151: return gzero;
1152: }
1153: int pari_compare_lg(GEN x, GEN y)
1154: {
1155: return lg(x)-lg(y);
1156: }
1157:
1158: /*
1.3 ! noro 1159: * Calcule la liste des sous groupes de (\ZZ/m\ZZ)^* d'ordre divisant p et
1.1 noro 1160: * retourne la liste de leurs elements
1161: */
1162: GEN
1163: listsousgroupes(long m, long p)
1164: {
1.3 ! noro 1165: gpmem_t ltop = avma;
! 1166: GEN zn, zns, lss, sg, res;
1.1 noro 1167: int k, card, i, phi;
1168: if (m == 2)
1169: {
1170: res = cgetg(2, t_VEC);
1171: sg = cgetg(2, t_VECSMALL);
1172: res[1] = (long) sg;
1173: sg[1] = 1;
1174: return res;
1175: }
1.3 ! noro 1176: zn = znstar(stoi(m));
! 1177: phi = itos((GEN) zn[1]);
! 1178: zns = znstar_small(zn);
! 1179: lss = subgrouplist((GEN) zn[2], NULL);
1.1 noro 1180: res = cgetg(lg(lss), t_VEC);
1181: for (k = 1, i = lg(lss) - 1; i >= 1; i--)
1182: {
1.3 ! noro 1183: gpmem_t av;
1.1 noro 1184: av = avma;
1185: card = phi / itos(det((GEN) lss[i]));
1186: avma = av;
1187: if (p % card == 0)
1.3 ! noro 1188: res[k++] = (long) znstar_hnf_elts(zns,(GEN)lss[i]);
1.1 noro 1189: }
1190: setlg(res,k);
1191: res=gen_sort(res,0,&pari_compare_lg);
1192: return gerepileupto(ltop,res);
1193: }
1194:
1195: GEN
1196: fixedfieldpolsigma(GEN sigma, GEN p, GEN Tp, GEN sym, GEN deg, long g)
1197: {
1.3 ! noro 1198: gpmem_t ltop=avma;
1.1 noro 1199: long i, j, npows;
1200: GEN s, f, pows;
1201: sigma=lift(gmul(sigma,gmodulsg(1,p)));
1202: f=polx[varn(sigma)];
1203: s=zeropol(varn(sigma));
1204: for(j=1;j<lg(sym);j++)
1205: if(sym[j])
1206: {
1207: s=FpX_add(s,FpX_Fp_mul(FpXQ_pow(f,stoi(deg[j]),Tp,p),stoi(sym[j]),p),p);
1208: }
1209: npows = brent_kung_optpow(lgef(Tp)-4,g-1);
1210: pows = FpXQ_powers(sigma,npows,Tp,p);
1211: for(i=2; i<=g;i++)
1212: {
1213: f=FpX_FpXQV_compo(f,pows,Tp,p);
1214: for(j=1;j<lg(sym);j++)
1215: if(sym[j])
1216: {
1217: s=FpX_add(s,FpX_Fp_mul(FpXQ_pow(f,stoi(deg[j]),Tp,p),stoi(sym[j]),p),p);
1218: }
1219: }
1220: return gerepileupto(ltop, s);
1221: }
1222:
1223: GEN
1224: fixedfieldfactmod(GEN Sp, GEN p, GEN Tmod)
1225: {
1226: long i;
1227: long l=lg(Tmod);
1228: GEN F=cgetg(l,t_VEC);
1229: for(i=1;i<l;i++)
1.3 ! noro 1230: F[i]=(long)FpXQ_minpoly(FpX_res(Sp,(GEN) Tmod[i],p), (GEN) Tmod[i],p);
1.1 noro 1231: return F;
1232: }
1233:
1234: GEN
1235: fixedfieldnewtonsumaut(GEN sigma, GEN p, GEN Tp, GEN e, long g)
1236: {
1.3 ! noro 1237: gpmem_t ltop=avma;
1.1 noro 1238: long i;
1239: GEN s,f,V;
1240: long rt;
1241: sigma=lift(gmul(sigma,gmodulsg(1,p)));
1242: f=polx[varn(sigma)];
1243: rt=brent_kung_optpow(lgef(Tp)-4,g-1);
1244: V=FpXQ_powers(sigma,rt,Tp,p);
1245: s=FpXQ_pow(f,e,Tp,p);
1246: for(i=2; i<=g;i++)
1247: {
1248: f=FpX_FpXQV_compo(f,V,Tp,p);
1249: s=FpX_add(s,FpXQ_pow(f,e,Tp,p),p);
1250: }
1251: return gerepileupto(ltop, s);
1252: }
1253:
1254: GEN
1255: fixedfieldnewtonsum(GEN O, GEN L, GEN mod, GEN e)
1256: {
1257: long f,g;
1258: long i,j;
1259: GEN PL;
1260: f=lg(O)-1;
1261: g=lg(O[1])-1;
1262: PL=cgetg(lg(O), t_COL);
1263: for(i=1; i<=f; i++)
1264: {
1.3 ! noro 1265: gpmem_t ltop=avma;
1.1 noro 1266: GEN s=gzero;
1267: for(j=1; j<=g; j++)
1268: s=addii(s,powmodulo((GEN)L[mael(O,i,j)],e,mod));
1269: PL[i]=lpileupto(ltop,modii(s,mod));
1270: }
1271: return PL;
1272: }
1273:
1274: GEN
1275: fixedfieldpol(GEN O, GEN L, GEN mod, GEN sym, GEN deg)
1276: {
1.3 ! noro 1277: gpmem_t ltop=avma;
1.1 noro 1278: long i;
1279: GEN S=gzero;
1280: for(i=1;i<lg(sym);i++)
1281: if (sym[i])
1282: S=gadd(S,gmulsg(sym[i],fixedfieldnewtonsum(O, L, mod, stoi(deg[i]))));
1283: S=FpV_red(S,mod);
1284: return gerepileupto(ltop, S);
1285: }
1286:
1287: static long
1288: fixedfieldtests(GEN LN, long n)
1289: {
1290: long i,j,k;
1291: for (i=1;i<lg(LN[1]);i++)
1292: for(j=i+1;j<lg(LN[1]);j++)
1293: {
1294: for(k=1;k<=n;k++)
1295: if (cmpii(gmael(LN,k,j),gmael(LN,k,i)))
1296: break;
1297: if (k>n)
1298: return 0;
1299: }
1300: return 1;
1301: }
1302:
1303: static long
1304: fixedfieldtest(GEN V)
1305: {
1306: long i,j;
1307: for (i=1;i<lg(V);i++)
1308: for(j=i+1;j<lg(V);j++)
1309: if (!cmpii((GEN)V[i],(GEN)V[j]))
1310: return 0;
1311: return 1;
1312: }
1313:
1314: void
1315: debug_surmer(char *s,GEN S, long n)
1316: {
1317: long l=lg(S);
1318: setlg(S,n+1);
1319: fprintferr(s,S);
1320: setlg(S,l);
1321: }
1322:
1323: /*We try hard to find a polynomial R squarefree mod p.
1324: Unfortunately, it may be too much asked.
1325: A less bugprone way would be to only ask it to be squarefree mod p^e
1326: with e not too big. Most of the code is here, but I lack the theoretical
1327: knowledge to make it to work smoothly.B.A.
1328: */
1329: static GEN
1330: fixedfieldsurmer(GEN O, GEN L, GEN mod, GEN l, GEN p, GEN S, GEN deg, long v, GEN LN,long n)
1331: {
1332: long i;
1333: GEN V;
1334: V=(GEN)LN[n];
1335: for (i=1;i<lg(S);i++)
1336: S[i]=0;
1337: S[n]=1;
1338: for (i=0;i<2*n-1;i++)
1339: {
1340: long k;
1341: if (fixedfieldtest(V))
1342: {
1.3 ! noro 1343: gpmem_t av=avma;
1.1 noro 1344: GEN s=fixedfieldpol(O,L,mod,S,deg);
1345: GEN P=FpV_roots_to_pol(s,mod,v);
1346: P=FpX_center(P,mod);
1347: if (p==gun || FpX_is_squarefree(P,p))
1348: {
1349: GEN V;
1350: if (DEBUGLEVEL>=4)
1351: debug_surmer("FixedField: Sym: %Z\n",S,n);
1352: V=cgetg(3,t_VEC);
1353: V[1]=lcopy(s);/*do not swap*/
1354: V[2]=lcopy(P);
1355: return V;
1356: }
1357: else
1358: {
1359: if (DEBUGLEVEL>=6)
1360: debug_surmer("FixedField: bad mod: %Z\n",S,n);
1361: avma=av;
1362: }
1363: }
1364: else
1365: {
1366: if (DEBUGLEVEL>=6)
1367: debug_surmer("FixedField: Tested: %Z\n",S,n);
1368: }
1369: k=1+(i%n);
1370: S[k]++;
1371: V=FpV_red(gadd(V,(GEN)LN[k]),l);
1372: }
1373: return NULL;
1374: }
1375:
1376: GEN
1377: fixedfieldsympol(GEN O, GEN L, GEN mod, GEN l, GEN p, GEN S, GEN deg, long v)
1378: {
1.3 ! noro 1379: gpmem_t ltop=avma;
1.1 noro 1380: GEN V=NULL;
1381: GEN LN=cgetg(lg(L),t_VEC);
1382: GEN Ll=FpV_red(L,l);
1383: long n,i;
1384: for(i=1;i<lg(deg);i++)
1385: deg[i]=0;
1386: for(n=1,i=1;i<lg(L);i++)
1387: {
1388: long j;
1389: LN[n]=(long)fixedfieldnewtonsum(O,Ll,l,stoi(i));
1390: for(j=2;j<lg(LN[n]);j++)
1391: if(cmpii(gmael(LN,n,j),gmael(LN,n,1)))
1392: break;
1393: if(j==lg(LN[n]) && j>2)
1394: continue;
1395: if (DEBUGLEVEL>=6) fprintferr("FixedField: LN[%d]=%Z\n",n,LN[n]);
1396: deg[n]=i;
1397: if (fixedfieldtests(LN,n))
1398: {
1399: V=fixedfieldsurmer(O,L,mod,l,p,S,deg,v,LN,n);
1400: if (V)
1401: break;
1402: }
1403: n++;
1404: }
1405: if (DEBUGLEVEL>=4)
1406: {
1407: long l=lg(deg);
1408: setlg(deg,n);
1409: fprintferr("FixedField: Computed degrees: %Z\n",deg);
1410: setlg(deg,l);
1411: }
1412: if (!V)
1413: err(talker, "prime too small in fixedfield");
1414: return gerepileupto(ltop,V);
1415: }
1416: /*
1417: * Calcule l'inclusion de R dans T i.e. S telque T|R\compo S
1418: * Ne recopie pas PL.
1419: */
1420: GEN
1421: fixedfieldinclusion(GEN O, GEN PL)
1422: {
1423: GEN S;
1424: int i, j;
1425: S = cgetg((lg(O) - 1) * (lg(O[1]) - 1) + 1, t_COL);
1426: for (i = 1; i < lg(O); i++)
1427: for (j = 1; j < lg(O[i]); j++)
1428: S[((GEN *) O)[i][j]] = PL[i];
1429: return S;
1430: }
1431:
1432: /*Usually mod is bigger than than den so there is no need to reduce it.*/
1433: GEN
1434: vandermondeinversemod(GEN L, GEN T, GEN den, GEN mod)
1435: {
1.3 ! noro 1436: gpmem_t av;
1.1 noro 1437: int i, j, n = lg(L);
1438: long x = varn(T);
1439: GEN M, P, Tp;
1440: M = cgetg(n, t_MAT);
1441: av=avma;
1442: Tp = gclone(FpX_red(deriv(T, x),mod)); /*clone*/
1443: avma=av;
1444: for (i = 1; i < n; i++)
1445: {
1446: GEN z;
1447: av = avma;
1448: z = mpinvmod(FpX_eval(Tp, (GEN) L[i],mod),mod);
1449: z = modii(mulii(den,z),mod);
1450: P = FpX_Fp_mul(FpX_div(T, deg1pol(gun,negi((GEN) L[i]),x),mod), z, mod);
1451: M[i] = lgetg(n, t_COL);
1452: for (j = 1; j < n; j++)
1453: mael(M,i,j) = lcopy((GEN) P[1 + j]);
1454: M[i] = lpileupto(av,(GEN)M[i]);
1455: }
1456: gunclone(Tp); /*unclone*/
1457: return M;
1458: }
1459: /* Calcule le polynome associe a un vecteur conjugue v*/
1460: static GEN
1461: vectopol(GEN v, GEN M, GEN den , GEN mod, long x)
1462: {
1.3 ! noro 1463: long n = lg(v), i, k;
! 1464: gpmem_t av;
1.1 noro 1465: GEN z = cgetg(n+1,t_POL),p1,p2,mod2;
1466: av=avma;
1467: mod2=gclone(shifti(mod,-1));/*clone*/
1468: avma=av;
1469: z[1]=evalsigne(1)|evalvarn(x)|evallgef(n+1);
1470: for (i=2; i<=n; i++)
1471: {
1472: p1=gzero; av=avma;
1473: for (k=1; k<n; k++)
1474: {
1475: p2=mulii(gcoeff(M,i-1,k),(GEN)v[k]);
1476: p1=addii(p1,p2);
1477: }
1478: p1=modii(p1,mod);
1479: if (cmpii(p1,mod2)>0) p1=subii(p1,mod);
1480: p1=gdiv(p1,den);
1481: z[i]=lpileupto(av,p1);
1482: }
1483: gunclone(mod2);/*unclone*/
1484: return normalizepol_i(z,n+1);
1485: }
1486: /* Calcule le polynome associe a une permutation des racines*/
1487: static GEN
1488: permtopol(GEN p, GEN L, GEN M, GEN den, GEN mod, long x)
1489: {
1.3 ! noro 1490: long n = lg(L), i, k;
! 1491: gpmem_t av;
1.1 noro 1492: GEN z = cgetg(n+1,t_POL),p1,p2,mod2;
1493: if (lg(p) != n) err(talker,"incorrect permutation in permtopol");
1494: av=avma;
1495: mod2=gclone(shifti(mod,-1)); /*clone*/
1496: avma=av;
1497: z[1]=evalsigne(1)|evalvarn(x)|evallgef(n+1);
1498: for (i=2; i<=n; i++)
1499: {
1500: p1=gzero; av=avma;
1501: for (k=1; k<n; k++)
1502: {
1503: p2=mulii(gcoeff(M,i-1,k),(GEN)L[p[k]]);
1504: p1=addii(p1,p2);
1505: }
1506: p1=modii(p1,mod);
1507: if (cmpii(p1,mod2)>0) p1=subii(p1,mod);
1508: p1=gdiv(p1,den);
1509: z[i]=lpileupto(av,p1);
1510: }
1511: gunclone(mod2); /*unclone*/
1512: return normalizepol_i(z,n+1);
1513: }
1514:
1515: GEN
1516: galoisgrouptopol( GEN res, GEN L, GEN M, GEN den, GEN mod, long v)
1517: {
1518: GEN aut = cgetg(lg(res), t_COL);
1519: long i;
1520: for (i = 1; i < lg(res); i++)
1521: {
1522: if (DEBUGLEVEL>=6)
1523: fprintferr("%d ",i);
1524: aut[i] = (long) permtopol((GEN) res[i], L, M, den, mod, v);
1525: }
1526: return aut;
1527: }
1528:
1529: /*
1530: * Casse l'orbite O en sous-orbite d'ordre premier correspondant a des
1531: * puissance de l'element
1532: */
1533: GEN
1534: splitorbite(GEN O)
1535: {
1.3 ! noro 1536: gpmem_t lbot, ltop = avma;
1.1 noro 1537: int i, n;
1.3 ! noro 1538: GEN fc, res;
1.1 noro 1539: n = lg(O[1]) - 1;
1.3 ! noro 1540: fc = decomp_primary_small(n);
1.1 noro 1541: lbot = avma;
1542: res = cgetg(3, t_VEC);
1543: res[1] = lgetg(lg(fc), t_VEC);
1544: res[2] = lgetg(lg(fc), t_VECSMALL);
1545: for (i = 1; i < lg(fc); i++)
1546: {
1.3 ! noro 1547: mael(res,1,lg(fc) - i) = (long)cyc_powtoperm(O, n / fc[i]);
! 1548: mael(res,2,lg(fc) - i) = fc[i];
1.1 noro 1549: }
1550: if ( DEBUGLEVEL>=4 )
1551: fprintferr("GaloisConj:splitorbite: %Z\n",res);
1552: return gerepile(ltop, lbot, res);
1553: }
1554:
1555: /*
1556: * structure qui contient l'analyse du groupe de Galois par etude des degres
1557: * de Frobenius:
1558: *
1559: * p: nombre premier a relever deg: degre du lift à effectuer pgp: plus grand
1560: * diviseur premier du degre de T.
1561: * l: un nombre premier tel que T est totalement decompose
1562: * modulo l ppp: plus petit diviseur premier du degre de T. primepointer:
1563: * permet de calculer les premiers suivant p.
1564: */
1.3 ! noro 1565: enum ga_code {ga_all_normal=1,ga_ext_2=2,ga_non_wss=4};
1.1 noro 1566: struct galois_analysis
1567: {
1568: long p;
1569: long deg;
1570: long ord;
1571: long l;
1572: long ppp;
1573: long p4;
1.3 ! noro 1574: enum ga_code group;
1.1 noro 1575: byteptr primepointer;
1576: };
1577: void
1578: galoisanalysis(GEN T, struct galois_analysis *ga, long calcul_l, long karma_type)
1579: {
1.3 ! noro 1580: gpmem_t ltop=avma;
1.1 noro 1581: long n,p;
1582: long i;
1.3 ! noro 1583: enum k_code {k_amoeba=0,k_snake=1,k_fish=2,k_bird=4,k_rodent=6,k_dog=8,k_human=9,k_cat=12} karma;
1.1 noro 1584: long group,linf;
1585: /*TODO: complete the table to at least 200*/
1586: const int prim_nonss_orders[]={36,48,56,60,72,75,80,96,108,120,132,0};
1587: GEN F,Fp,Fe,Fpe,O;
1.3 ! noro 1588: long min_prime,np;
1.1 noro 1589: long order,phi_order;
1590: long plift,nbmax,nbtest,deg;
1591: byteptr primepointer,pp;
1.3 ! noro 1592:
! 1593: if (!ZX_is_squarefree(T))
! 1594: err(talker, "Polynomial not squarefree in galoisinit");
! 1595: if (DEBUGLEVEL >= 1) (void)timer2();
1.1 noro 1596: n = degpol(T);
1597: if (!karma_type) karma_type=n;
1598: else err(warner,"entering black magic computation");
1599: O = cgetg(n+1,t_VECSMALL);
1600: for(i=1;i<=n;i++) O[i]=0;
1601: F = factor(stoi(n));
1602: Fp=gtovecsmall((GEN)F[1]);
1603: Fe=gtovecsmall((GEN)F[2]);
1604: np=lg(Fp)-1;
1605: Fpe=cgetg(np+1, t_VECSMALL);
1606: for (i = 1; i < lg(Fpe); i++)
1607: Fpe[i] = itos(powgi(gmael(F,1,i), gmael(F,2,i)));
1608: /*In this part, we study the cardinal of the group to have an information
1609: about the orders, so if we are unlucky we can continue.*/
1610:
1611: /*Are there non WSS groups of this order ?*/
1612: group=0;
1613: for(i=0;prim_nonss_orders[i];i++)
1614: if (n%prim_nonss_orders[i] == 0)
1615: {
1616: group |= ga_non_wss;
1617: break;
1618: }
1619: if ( n>12 && n%12 == 0 )
1620: {
1621: /*We need to know the greatest prime dividing n/12*/
1622: if ( Fp[np] == 3 && Fe[np] == 1 )
1623: group |= ga_ext_2;
1624: }
1625: phi_order = 1;
1626: order = 1;
1627: for (i = np; i > 0; i--)
1628: {
1629: p = Fp[i];
1630: if (phi_order % p != 0)
1631: {
1632: order *= p;
1633: phi_order *= p - 1;
1634: }
1635: else
1636: {
1637: group |= ga_all_normal;
1638: break;
1639: }
1640: if (Fe[i]>1)
1641: break;
1642: }
1643: /*Now, we study the orders of the Frobenius elements*/
1.3 ! noro 1644: min_prime=n*max((long)(BITS_IN_LONG-bfffo(n)-4),2);
1.1 noro 1645: plift = 0;
1646: nbmax = 8+(n>>1);
1647: nbtest = 0; karma = k_amoeba;
1648: deg = Fp[np];
1649: for (p = 0, pp = primepointer = diffptr;
1650: (plift == 0
1651: || (nbtest < nbmax && (nbtest <=8 || order < (n>>1)))
1652: || (n == 24 && O[6] == 0 && O[4] == 0)
1653: || ((group&ga_non_wss) && order == Fp[np]))
1654: && (nbtest < 3 * nbmax || !(group&ga_non_wss)) ;)
1655: {
1.3 ! noro 1656: gpmem_t av;
1.1 noro 1657: GEN ip,FS,p1;
1658: long o,norm_o=1;
1.3 ! noro 1659:
! 1660: NEXT_PRIME_VIADIFF_CHECK(p,primepointer);
1.1 noro 1661: /*discard small primes*/
1.3 ! noro 1662: if (p <= min_prime)
1.1 noro 1663: continue;
1664: ip=stoi(p);
1665: if (!FpX_is_squarefree(T,ip))
1666: continue;
1667: nbtest++;
1668: av=avma;
1669: FS=(GEN)simplefactmod(T,ip)[1];
1670: p1=(GEN)FS[1];
1671: for(i=2;i<lg(FS);i++)
1672: if (cmpii(p1,(GEN)FS[i]))
1673: break;
1674: if (i<lg(FS))
1675: {
1676: avma = ltop;
1677: if (DEBUGLEVEL >= 2)
1678: fprintferr("GaloisAnalysis:non Galois for p=%ld\n", p);
1679: ga->p = p;
1680: ga->deg = 0;
1681: return; /* Not a Galois polynomial */
1682: }
1683: o=n/(lg(FS)-1);
1684: avma=av;
1685: if (!O[o])
1686: O[o]=p;
1687: if (o % deg == 0)
1688: {
1689: /*We try to find a power of the Frobenius which generate
1690: a normal subgroup just by looking at the order.*/
1691: if (o * Fp[1] >= n)
1692: /*Subgroup of smallest index are normal*/
1693: norm_o = o;
1694: else
1695: {
1696: norm_o = 1;
1697: for (i = np; i > 0; i--)
1698: {
1699: if (o % Fpe[i] == 0)
1700: norm_o *= Fpe[i];
1701: else
1702: break;
1703: }
1704: }
1705: if (norm_o != 1)
1706: {
1707: if (!(group&ga_all_normal) || o > order ||
1708: (o == order && (plift == 0 || norm_o > deg
1709: || (norm_o == deg && cgcd(p-1,karma_type)>karma ))))
1710: {
1711: deg = norm_o;
1712: order = o;
1713: plift = p;
1.3 ! noro 1714: karma=(enum k_code)cgcd(p-1,karma_type);
1.1 noro 1715: pp = primepointer;
1716: group |= ga_all_normal;
1717: }
1718: }
1719: else if (!(group&ga_all_normal) && (plift == 0 || o > order
1720: || ( o == order && cgcd(p-1,karma_type)>karma )))
1721: {
1722: order = o;
1723: plift = p;
1.3 ! noro 1724: karma=(enum k_code)cgcd(p-1,karma_type);
1.1 noro 1725: pp = primepointer;
1726: }
1727: }
1728: if (DEBUGLEVEL >= 5)
1729: fprintferr("GaloisAnalysis:Nbtest=%ld,p=%ld,o=%ld,n_o=%d,best p=%ld,ord=%ld,k=%ld\n",
1730: nbtest, p, o, norm_o, plift, order,karma);
1731: }
1732: /* This is to avoid looping on non-wss group.
1733: To be checked for large groups. */
1734: /*Would it be better to disable this check if we are in a good case ?
1735: * (ga_all_normal and !(ga_ext_2) (e.g. 60)) ?*/
1736: if (plift == 0 || ((group&ga_non_wss) && order == Fp[np]))
1737: {
1738: deg = 0;
1739: err(warner, "Galois group almost certainly not weakly super solvable");
1740: }
1741: /*linf=(n*(n-1))>>2;*/
1742: linf=n;
1743: if (calcul_l && O[1]<=linf)
1744: {
1.3 ! noro 1745: gpmem_t av;
1.1 noro 1746: long l=0;
1747: /*we need a totally splited prime l*/
1748: av = avma;
1749: while (l == 0)
1750: {
1751: long nb;
1.3 ! noro 1752:
! 1753: NEXT_PRIME_VIADIFF_CHECK(p,primepointer);
1.1 noro 1754: if (p<=linf) continue;
1755: nb=FpX_nbroots(T,stoi(p));
1756: if (nb == n)
1757: l = p;
1758: else if (nb && FpX_is_squarefree(T,stoi(p)))
1759: {
1760: avma = ltop;
1761: if (DEBUGLEVEL >= 2)
1762: fprintferr("GaloisAnalysis:non Galois for p=%ld\n", p);
1763: ga->p = p;
1764: ga->deg = 0;
1765: return; /* Not a Galois polynomial */
1766: }
1767: avma = av;
1768: }
1769: O[1]=l;
1770: }
1771: ga->p = plift;
1.3 ! noro 1772: ga->group = (enum ga_code)group;
1.1 noro 1773: ga->deg = deg;
1774: ga->ord = order;
1775: ga->l = O[1];
1776: ga->primepointer = pp;
1777: ga->ppp = Fp[1];
1778: ga->p4 = O[4];
1779: if (DEBUGLEVEL >= 4)
1780: fprintferr("GaloisAnalysis:p=%ld l=%ld group=%ld deg=%ld ord=%ld\n",
1781: plift, O[1], group, deg, order);
1782: if (DEBUGLEVEL >= 1)
1783: msgtimer("galoisanalysis()");
1784: avma = ltop;
1785: }
1786:
1787: /* Groupe A4 */
1788: GEN
1789: a4galoisgen(GEN T, struct galois_test *td)
1790: {
1.3 ! noro 1791: gpmem_t ltop = avma, av, av2;
1.1 noro 1792: long i, j, k;
1793: long n;
1794: long N, hop = 0;
1795: long **O;
1796: GEN *ar, **mt; /* tired of casting */
1797: GEN t, u;
1798: GEN res, orb, ry;
1799: GEN pft, pfu, pfv;
1800: n = degpol(T);
1801: res = cgetg(3, t_VEC);
1802: ry = cgetg(4, t_VEC);
1803: res[1] = (long) ry;
1804: pft = cgetg(n + 1, t_VECSMALL);
1805: pfu = cgetg(n + 1, t_VECSMALL);
1806: pfv = cgetg(n + 1, t_VECSMALL);
1807: ry[1] = (long) pft;
1808: ry[2] = (long) pfu;
1809: ry[3] = (long) pfv;
1810: ry = cgetg(4, t_VECSMALL);
1811: ry[1] = 2;
1812: ry[2] = 2;
1813: ry[3] = 3;
1814: res[2] = (long) ry;
1815: av = avma;
1816: ar = (GEN *) alloue_ardoise(n, 1 + lg(td->ladic));
1817: mt = (GEN **) td->PV[td->ordre[n]];
1818: t = cgetg(n + 1, t_VECSMALL) + 1; /* Sorry for this hack */
1819: u = cgetg(n + 1, t_VECSMALL) + 1; /* too lazy to correct */
1820: av2 = avma;
1821: N = itos(gdiv(mpfact(n), mpfact(n >> 1))) >> (n >> 1);
1822: if (DEBUGLEVEL >= 4)
1823: fprintferr("A4GaloisConj:I will test %ld permutations\n", N);
1824: avma = av2;
1825: for (i = 0; i < n; i++)
1826: t[i] = i + 1;
1827: for (i = 0; i < N; i++)
1828: {
1829: GEN g;
1830: int a, x, y;
1831: if (i == 0)
1832: {
1833: gaffect(gzero, ar[(n - 2) >> 1]);
1834: for (k = n - 2; k > 2; k -= 2)
1835: addiiz(ar[k >> 1], addii(mt[k + 1][k + 2], mt[k + 2][k + 1]),
1836: ar[(k >> 1) - 1]);
1837: }
1838: else
1839: {
1840: x = i;
1841: y = 1;
1842: do
1843: {
1844: y += 2;
1845: a = x%y;
1846: x = x/y;
1847: }
1848: while (!a);
1849: switch (y)
1850: {
1851: case 3:
1852: x = t[2];
1853: if (a == 1)
1854: {
1855: t[2] = t[1];
1856: t[1] = x;
1857: }
1858: else
1859: {
1860: t[2] = t[0];
1861: t[0] = x;
1862: }
1863: break;
1864: case 5:
1865: x = t[0];
1866: t[0] = t[2];
1867: t[2] = t[1];
1868: t[1] = x;
1869: x = t[4];
1870: t[4] = t[4 - a];
1871: t[4 - a] = x;
1872: addiiz(ar[2], addii(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
1873: break;
1874: case 7:
1875: x = t[0];
1876: t[0] = t[4];
1877: t[4] = t[3];
1878: t[3] = t[1];
1879: t[1] = t[2];
1880: t[2] = x;
1881: x = t[6];
1882: t[6] = t[6 - a];
1883: t[6 - a] = x;
1884: addiiz(ar[3], addii(mt[t[6]][t[7]], mt[t[7]][t[6]]), ar[2]);
1885: addiiz(ar[2], addii(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
1886: break;
1887: case 9:
1888: x = t[0];
1889: t[0] = t[6];
1890: t[6] = t[5];
1891: t[5] = t[3];
1892: t[3] = x;
1893: x = t[4];
1894: t[4] = t[1];
1895: t[1] = x;
1896: x = t[8];
1897: t[8] = t[8 - a];
1898: t[8 - a] = x;
1899: addiiz(ar[4], addii(mt[t[8]][t[9]], mt[t[9]][t[8]]), ar[3]);
1900: addiiz(ar[3], addii(mt[t[6]][t[7]], mt[t[7]][t[6]]), ar[2]);
1901: addiiz(ar[2], addii(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
1902: break;
1903: default:
1904: y--;
1905: x = t[0];
1906: t[0] = t[2];
1907: t[2] = t[1];
1908: t[1] = x;
1909: for (k = 4; k < y; k += 2)
1910: {
1911: int j;
1912: x = t[k];
1913: for (j = k; j > 0; j--)
1914: t[j] = t[j - 1];
1915: t[0] = x;
1916: }
1917: x = t[y];
1918: t[y] = t[y - a];
1919: t[y - a] = x;
1920: for (k = y; k > 2; k -= 2)
1921: addiiz(ar[k >> 1],
1922: addii(mt[t[k]][t[k + 1]], mt[t[k + 1]][t[k]]),
1923: ar[(k >> 1) - 1]);
1924: }
1925: }
1926: g = addii(ar[1], addii(addii(mt[t[0]][t[1]], mt[t[1]][t[0]]),
1927: addii(mt[t[2]][t[3]], mt[t[3]][t[2]])));
1928: if (padicisint(g, td))
1929: {
1930: for (k = 0; k < n; k += 2)
1931: {
1932: pft[t[k]] = t[k + 1];
1933: pft[t[k + 1]] = t[k];
1934: }
1935: if (verifietest(pft, td))
1936: break;
1937: else
1938: hop++;
1939: }
1940: avma = av2;
1941: }
1942: if (i == N)
1943: {
1944: avma = ltop;
1945: if (DEBUGLEVEL >= 1 && hop)
1946: fprintferr("A4GaloisConj: %ld hop sur %ld iterations\n", hop, N);
1947: return gzero;
1948: }
1949: if (DEBUGLEVEL >= 1 && hop)
1950: fprintferr("A4GaloisConj: %ld hop sur %ld iterations\n", hop, N);
1951: N = itos(gdiv(mpfact(n >> 1), mpfact(n >> 2))) >> 1;
1952: avma = av2;
1953: if (DEBUGLEVEL >= 4)
1954: fprintferr("A4GaloisConj:sigma=%Z \n", pft);
1955: for (i = 0; i < N; i++)
1956: {
1957: GEN g;
1958: int a, x, y;
1959: if (i == 0)
1960: {
1961: for (k = 0; k < n; k += 4)
1962: {
1963: u[k + 3] = t[k + 3];
1964: u[k + 2] = t[k + 1];
1965: u[k + 1] = t[k + 2];
1966: u[k] = t[k];
1967: }
1968: }
1969: else
1970: {
1971: x = i;
1972: y = -2;
1973: do
1974: {
1975: y += 4;
1976: a = x%y;
1977: x = x/y;
1978: }
1979: while (!a);
1980: x = u[2];
1981: u[2] = u[0];
1982: u[0] = x;
1983: switch (y)
1984: {
1985: case 2:
1986: break;
1987: case 6:
1988: x = u[4];
1989: u[4] = u[6];
1990: u[6] = x;
1991: if (!(a & 1))
1992: {
1993: a = 4 - (a >> 1);
1994: x = u[6];
1995: u[6] = u[a];
1996: u[a] = x;
1997: x = u[4];
1998: u[4] = u[a - 2];
1999: u[a - 2] = x;
2000: }
2001: break;
2002: case 10:
2003: x = u[6];
2004: u[6] = u[3];
2005: u[3] = u[2];
2006: u[2] = u[4];
2007: u[4] = u[1];
2008: u[1] = u[0];
2009: u[0] = x;
2010: if (a >= 3)
2011: a += 2;
2012: a = 8 - a;
2013: x = u[10];
2014: u[10] = u[a];
2015: u[a] = x;
2016: x = u[8];
2017: u[8] = u[a - 2];
2018: u[a - 2] = x;
2019: break;
2020: }
2021: }
2022: g = gzero;
2023: for (k = 0; k < n; k += 2)
2024: g = addii(g, addii(mt[u[k]][u[k + 1]], mt[u[k + 1]][u[k]]));
2025: if (padicisint(g, td))
2026: {
2027: for (k = 0; k < n; k += 2)
2028: {
2029: pfu[u[k]] = u[k + 1];
2030: pfu[u[k + 1]] = u[k];
2031: }
2032: if (verifietest(pfu, td))
2033: break;
2034: else
2035: hop++;
2036: }
2037: avma = av2;
2038: }
2039: if (i == N)
2040: {
2041: avma = ltop;
2042: return gzero;
2043: }
2044: if (DEBUGLEVEL >= 1 && hop)
2045: fprintferr("A4GaloisConj: %ld hop sur %ld iterations\n", hop, N);
2046: if (DEBUGLEVEL >= 4)
2047: fprintferr("A4GaloisConj:tau=%Z \n", pfu);
2048: avma = av2;
2049: orb = cgetg(3, t_VEC);
2050: orb[1] = (long) pft;
2051: orb[2] = (long) pfu;
2052: if (DEBUGLEVEL >= 4)
2053: fprintferr("A4GaloisConj:orb=%Z \n", orb);
1.3 ! noro 2054: O = (long**) vecperm_orbits(orb, 12);
1.1 noro 2055: if (DEBUGLEVEL >= 4)
2056: fprintferr("A4GaloisConj:O=%Z \n", O);
2057: av2 = avma;
2058: for (j = 0; j < 2; j++)
2059: {
2060: pfv[O[1][1]] = O[2][1];
2061: pfv[O[1][2]] = O[2][3 + j];
2062: pfv[O[1][3]] = O[2][4 - (j << 1)];
2063: pfv[O[1][4]] = O[2][2 + j];
2064: for (i = 0; i < 4; i++)
2065: {
2066: long x;
2067: GEN g;
2068: switch (i)
2069: {
2070: case 0:
2071: break;
2072: case 1:
2073: x = O[3][1];
2074: O[3][1] = O[3][2];
2075: O[3][2] = x;
2076: x = O[3][3];
2077: O[3][3] = O[3][4];
2078: O[3][4] = x;
2079: break;
2080: case 2:
2081: x = O[3][1];
2082: O[3][1] = O[3][4];
2083: O[3][4] = x;
2084: x = O[3][2];
2085: O[3][2] = O[3][3];
2086: O[3][3] = x;
2087: break;
2088: case 3:
2089: x = O[3][1];
2090: O[3][1] = O[3][2];
2091: O[3][2] = x;
2092: x = O[3][3];
2093: O[3][3] = O[3][4];
2094: O[3][4] = x;
2095: }
2096: pfv[O[2][1]] = O[3][1];
2097: pfv[O[2][3 + j]] = O[3][4 - j];
2098: pfv[O[2][4 - (j << 1)]] = O[3][2 + (j << 1)];
2099: pfv[O[2][2 + j]] = O[3][3 - j];
2100: pfv[O[3][1]] = O[1][1];
2101: pfv[O[3][4 - j]] = O[1][2];
2102: pfv[O[3][2 + (j << 1)]] = O[1][3];
2103: pfv[O[3][3 - j]] = O[1][4];
2104: g = gzero;
2105: for (k = 1; k <= n; k++)
2106: g = addii(g, mt[k][pfv[k]]);
2107: if (padicisint(g, td) && verifietest(pfv, td))
2108: {
2109: avma = av;
2110: if (DEBUGLEVEL >= 1)
2111: fprintferr("A4GaloisConj:%ld hop sur %d iterations max\n",
2112: hop, 10395 + 68);
2113: return res;
2114: }
2115: else
2116: hop++;
2117: avma = av2;
2118: }
2119: }
2120: /* Echec? */
2121: avma = ltop;
2122: return gzero;
2123: }
2124:
2125: /* Groupe S4 */
2126: static void
2127: s4makelift(GEN u, struct galois_lift *gl, GEN liftpow)
2128: {
2129: int i;
2130: liftpow[1] = (long) automorphismlift(u, gl, NULL);
2131: for (i = 2; i < lg(liftpow); i++)
2132: liftpow[i] = (long) FpXQ_mul((GEN) liftpow[i - 1], (GEN) liftpow[1],gl->TQ,gl->Q);
2133: }
2134: static long
2135: s4test(GEN u, GEN liftpow, struct galois_lift *gl, GEN phi)
2136: {
1.3 ! noro 2137: gpmem_t ltop=avma;
1.1 noro 2138: GEN res;
2139: int bl,i,d = lgef(u)-2;
1.3 ! noro 2140: if (DEBUGLEVEL >= 6) (void)timer2();
1.1 noro 2141: if ( !d ) return 0;
2142: res=(GEN)u[2];
2143: for (i = 1; i < d; i++)
2144: {
2145: if (lgef(liftpow[i])>2)
2146: res=addii(res,mulii(gmael(liftpow,i,2), (GEN) u[i + 2]));
2147: }
2148: res=modii(mulii(res,gl->den),gl->Q);
2149: if (cmpii(res,gl->gb->bornesol)>0
2150: && cmpii(res,subii(gl->Q,gl->gb->bornesol))<0)
2151: {
2152: avma=ltop;
2153: return 0;
2154: }
2155: res = (GEN) scalarpol((GEN)u[2],varn(u));
2156: for (i = 1; i < d ; i++)
2157: {
2158: GEN z;
2159: z = FpX_Fp_mul((GEN) liftpow[i], (GEN) u[i + 2],NULL);
2160: res = FpX_add(res,z ,gl->Q);
2161: }
2162: res = FpX_center(FpX_Fp_mul(res,gl->den,gl->Q), gl->Q);
2163: if (DEBUGLEVEL >= 6)
2164: msgtimer("s4test()");
2165: bl = poltopermtest(res, gl, phi);
2166: avma=ltop;
2167: return bl;
2168: }
2169: static GEN
2170: s4releveauto(GEN misom,GEN Tmod,GEN Tp,GEN p,long a1,long a2,long a3,long a4,long a5,long a6)
2171: {
1.3 ! noro 2172: gpmem_t ltop=avma;
1.1 noro 2173: GEN u1,u2,u3,u4,u5;
2174: GEN pu1,pu2,pu3,pu4;
2175: pu1=FpX_mul( (GEN) Tmod[a2], (GEN) Tmod[a1],p);
2176: u1 = FpX_chinese_coprime(gmael(misom,a1,a2),gmael(misom,a2,a1),
2177: (GEN) Tmod[a2], (GEN) Tmod[a1],pu1,p);
2178: pu2=FpX_mul( (GEN) Tmod[a4], (GEN) Tmod[a3],p);
2179: u2 = FpX_chinese_coprime(gmael(misom,a3,a4),gmael(misom,a4,a3),
2180: (GEN) Tmod[a4], (GEN) Tmod[a3],pu2,p);
2181: pu3=FpX_mul( (GEN) Tmod[a6], (GEN) Tmod[a5],p);
2182: u3 = FpX_chinese_coprime(gmael(misom,a5,a6),gmael(misom,a6,a5),
2183: (GEN) Tmod[a6], (GEN) Tmod[a5],pu3,p);
2184: pu4=FpX_mul(pu1,pu2,p);
2185: u4 = FpX_chinese_coprime(u1,u2,pu1,pu2,pu4,p);
2186: u5 = FpX_chinese_coprime(u4,u3,pu4,pu3,Tp,p);
2187: return gerepileupto(ltop,u5);
2188: }
2189: GEN
2190: s4galoisgen(struct galois_lift *gl)
2191: {
2192: struct galois_testlift gt;
1.3 ! noro 2193: gpmem_t av, ltop2, ltop = avma;
1.1 noro 2194: GEN Tmod, isom, isominv, misom;
2195: int i, j;
2196: GEN sg;
2197: GEN sigma, tau, phi;
2198: GEN res, ry;
2199: GEN pj;
2200: GEN p,Q,TQ,Tp;
2201: GEN bezoutcoeff, pauto, liftpow, aut;
2202:
2203: p = gl->p;
2204: Q = gl->Q;
2205: res = cgetg(3, t_VEC);
2206: ry = cgetg(5, t_VEC);
2207: res[1] = (long) ry;
2208: for (i = 1; i < lg(ry); i++)
2209: ry[i] = lgetg(lg(gl->L), t_VECSMALL);
2210: ry = cgetg(5, t_VECSMALL);
2211: res[2] = (long) ry;
2212: ry[1] = 2;
2213: ry[2] = 2;
2214: ry[3] = 3;
2215: ry[4] = 2;
2216: ltop2 = avma;
2217: sg = cgetg(7, t_VECSMALL);
2218: pj = cgetg(7, t_VECSMALL);
2219: sigma = cgetg(lg(gl->L), t_VECSMALL);
2220: tau = cgetg(lg(gl->L), t_VECSMALL);
2221: phi = cgetg(lg(gl->L), t_VECSMALL);
2222: for (i = 1; i < lg(sg); i++)
2223: sg[i] = i;
2224: Tp = FpX_red(gl->T,p);
2225: TQ = gl->TQ;
2226: Tmod = lift((GEN) factmod(gl->T, p)[1]);
2227: isom = cgetg(lg(Tmod), t_VEC);
2228: isominv = cgetg(lg(Tmod), t_VEC);
2229: misom = cgetg(lg(Tmod), t_MAT);
2230: aut=galoisdolift(gl, NULL);
2231: inittestlift(aut,Tmod, gl, >);
2232: bezoutcoeff = gt.bezoutcoeff;
2233: pauto = gt.pauto;
2234: for (i = 1; i < lg(pj); i++)
2235: pj[i] = 0;
2236: for (i = 1; i < lg(isom); i++)
2237: {
2238: misom[i] = lgetg(lg(Tmod), t_COL);
2239: isom[i] = (long) Fp_isom((GEN) Tmod[1], (GEN) Tmod[i], p);
2240: if (DEBUGLEVEL >= 6)
2241: fprintferr("S4GaloisConj:Computing isomorphisms %d:%Z\n", i,
2242: (GEN) isom[i]);
2243: isominv[i] = (long) Fp_inv_isom((GEN) isom[i], (GEN) Tmod[i],p);
2244: }
2245: for (i = 1; i < lg(isom); i++)
2246: for (j = 1; j < lg(isom); j++)
2247: mael(misom,i,j) = (long) FpX_FpXQ_compo((GEN) isominv[i],(GEN) isom[j],
2248: (GEN) Tmod[j],p);
2249: liftpow = cgetg(24, t_VEC);
2250: av = avma;
2251: for (i = 0; i < 3; i++)
2252: {
1.3 ! noro 2253: gpmem_t av2, avm1, avm2;
1.1 noro 2254: GEN u;
2255: int j1, j2, j3;
2256: GEN u1, u2, u3;
2257: if (i)
2258: {
2259: int x;
2260: x = sg[3];
2261: if (i == 1)
2262: {
2263: sg[3] = sg[2];
2264: sg[2] = x;
2265: }
2266: else
2267: {
2268: sg[3] = sg[1];
2269: sg[1] = x;
2270: }
2271: }
2272: u=s4releveauto(misom,Tmod,Tp,p,sg[1],sg[2],sg[3],sg[4],sg[5],sg[6]);
2273: s4makelift(u, gl, liftpow);
2274: av2 = avma;
2275: for (j1 = 0; j1 < 4; j1++)
2276: {
2277: u1 = FpX_add(FpXQ_mul((GEN) bezoutcoeff[sg[5]],
2278: (GEN) pauto[1 + j1],TQ,Q),
2279: FpXQ_mul((GEN) bezoutcoeff[sg[6]],
2280: (GEN) pauto[((-j1) & 3) + 1],TQ,Q),Q);
2281: avm1 = avma;
2282: for (j2 = 0; j2 < 4; j2++)
2283: {
2284: u2 = FpX_add(u1, FpXQ_mul((GEN) bezoutcoeff[sg[3]],
2285: (GEN) pauto[1 + j2],TQ,Q),NULL);
2286: u2 = FpX_add(u2, FpXQ_mul((GEN) bezoutcoeff[sg[4]], (GEN)
2287: pauto[((-j2) & 3) + 1],TQ,Q),Q);
2288: avm2 = avma;
2289: for (j3 = 0; j3 < 4; j3++)
2290: {
2291: u3 = FpX_add(u2, FpXQ_mul((GEN) bezoutcoeff[sg[1]],
2292: (GEN) pauto[1 + j3],TQ,Q),NULL);
2293: u3 = FpX_add(u3, FpXQ_mul((GEN) bezoutcoeff[sg[2]], (GEN)
2294: pauto[((-j3) & 3) + 1],TQ,Q),Q);
2295: if (DEBUGLEVEL >= 4)
2296: fprintferr("S4GaloisConj:Testing %d/3:%d/4:%d/4:%d/4:%Z\n",
2297: i, j1,j2, j3, sg);
2298: if (s4test(u3, liftpow, gl, sigma))
2299: {
2300: pj[1] = j3;
2301: pj[2] = j2;
2302: pj[3] = j1;
2303: goto suites4;
2304: }
2305: avma = avm2;
2306: }
2307: avma = avm1;
2308: }
2309: avma = av2;
2310: }
2311: avma = av;
2312: }
2313: avma = ltop;
2314: return gzero;
2315: suites4:
2316: if (DEBUGLEVEL >= 4)
2317: fprintferr("S4GaloisConj:sigma=%Z\n", sigma);
2318: if (DEBUGLEVEL >= 4)
2319: fprintferr("S4GaloisConj:pj=%Z\n", pj);
2320: avma = av;
2321: for (j = 1; j <= 3; j++)
2322: {
1.3 ! noro 2323: gpmem_t av2;
1.1 noro 2324: GEN u;
2325: int w, l;
2326: int z;
2327: z = sg[1]; sg[1] = sg[3]; sg[3] = sg[5]; sg[5] = z;
2328: z = sg[2]; sg[2] = sg[4]; sg[4] = sg[6]; sg[6] = z;
2329: z = pj[1]; pj[1] = pj[2]; pj[2] = pj[3]; pj[3] = z;
2330: for (l = 0; l < 2; l++)
2331: {
2332: u=s4releveauto(misom,Tmod,Tp,p,sg[1],sg[3],sg[2],sg[4],sg[5],sg[6]);
2333: s4makelift(u, gl, liftpow);
2334: av2 = avma;
2335: for (w = 0; w < 4; w += 2)
2336: {
1.3 ! noro 2337: gpmem_t av3;
1.1 noro 2338: GEN uu;
2339: pj[6] = (w + pj[3]) & 3;
2340: uu =FpX_add(FpXQ_mul((GEN) bezoutcoeff[sg[5]],
2341: (GEN) pauto[(pj[6] & 3) + 1],TQ,Q),
2342: FpXQ_mul((GEN) bezoutcoeff[sg[6]],
2343: (GEN) pauto[((-pj[6]) & 3) + 1],TQ,Q),Q);
2344: av3 = avma;
2345: for (i = 0; i < 4; i++)
2346: {
2347: GEN u;
2348: pj[4] = i;
2349: pj[5] = (i + pj[2] - pj[1]) & 3;
2350: if (DEBUGLEVEL >= 4)
2351: fprintferr("S4GaloisConj:Testing %d/3:%d/2:%d/2:%d/4:%Z:%Z\n",
2352: j - 1, w >> 1, l, i, sg, pj);
2353: u = FpX_add(uu, FpXQ_mul((GEN) pauto[(pj[4] & 3) + 1],
2354: (GEN) bezoutcoeff[sg[1]],TQ,Q),Q);
2355: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-pj[4]) & 3) + 1],
2356: (GEN) bezoutcoeff[sg[3]],TQ,Q),Q);
2357: u = FpX_add(u, FpXQ_mul((GEN) pauto[(pj[5] & 3) + 1],
2358: (GEN) bezoutcoeff[sg[2]],TQ,Q),Q);
2359: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-pj[5]) & 3) + 1],
2360: (GEN) bezoutcoeff[sg[4]],TQ,Q),Q);
2361: if (s4test(u, liftpow, gl, tau))
2362: goto suites4_2;
2363: avma = av3;
2364: }
2365: avma = av2;
2366: }
2367: z = sg[4];
2368: sg[4] = sg[3];
2369: sg[3] = z;
2370: pj[2] = (-pj[2]) & 3;
2371: avma = av;
2372: }
2373: }
2374: avma = ltop;
2375: return gzero;
2376: suites4_2:
2377: avma = av;
2378: {
2379: int abc, abcdef;
2380: GEN u;
1.3 ! noro 2381: gpmem_t av2;
1.1 noro 2382: abc = (pj[1] + pj[2] + pj[3]) & 3;
2383: abcdef = (((abc + pj[4] + pj[5] - pj[6]) & 3) >> 1);
2384: u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[4],sg[2],sg[5],sg[3],sg[6]);
2385: s4makelift(u, gl, liftpow);
2386: av2 = avma;
2387: for (j = 0; j < 8; j++)
2388: {
2389: int h, g, i;
2390: h = j & 3;
2391: g = abcdef + ((j & 4) >> 1);
2392: i = h + abc - g;
2393: u = FpXQ_mul((GEN) pauto[(g & 3) + 1],
2394: (GEN) bezoutcoeff[sg[1]],TQ,Q);
2395: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-g) & 3) + 1],
2396: (GEN) bezoutcoeff[sg[4]],TQ,Q),NULL);
2397: u = FpX_add(u, FpXQ_mul((GEN) pauto[(h & 3) + 1],
2398: (GEN) bezoutcoeff[sg[2]],TQ,Q),NULL);
2399: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-h) & 3) + 1],
2400: (GEN) bezoutcoeff[sg[5]],TQ,Q),NULL);
2401: u = FpX_add(u, FpXQ_mul((GEN) pauto[(i & 3) + 1],
2402: (GEN) bezoutcoeff[sg[3]],TQ,Q),NULL);
2403: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-i) & 3) + 1],
2404: (GEN) bezoutcoeff[sg[6]],TQ,Q),Q);
2405: if (DEBUGLEVEL >= 4)
2406: fprintferr("S4GaloisConj:Testing %d/8 %d:%d:%d\n",
2407: j, g & 3, h & 3, i & 3);
2408: if (s4test(u, liftpow, gl, phi))
2409: break;
2410: avma = av2;
2411: }
2412: }
2413: if (j == 8)
2414: {
2415: avma = ltop;
2416: return gzero;
2417: }
2418: for (i = 1; i < lg(gl->L); i++)
2419: {
2420: ((GEN **) res)[1][1][i] = sigma[tau[i]];
2421: ((GEN **) res)[1][2][i] = phi[sigma[tau[phi[i]]]];
2422: ((GEN **) res)[1][3][i] = phi[sigma[i]];
2423: ((GEN **) res)[1][4][i] = sigma[i];
2424: }
2425: avma = ltop2;
2426: return res;
2427: }
2428: struct galois_frobenius
2429: {
2430: long p;
2431: long fp;
2432: long deg;
2433: GEN Tmod;
2434: GEN psi;
2435: };
2436:
2437: /*Warning : the output of this function is not gerepileupto
2438: * compatible...*/
2439: static GEN
2440: galoisfindgroups(GEN lo, GEN sg, long f)
2441: {
1.3 ! noro 2442: gpmem_t ltop=avma;
1.1 noro 2443: GEN V,W;
2444: long i,j,k;
2445: V=cgetg(lg(lo),t_VEC);
2446: for(j=1,i=1;i<lg(lo);i++)
2447: {
1.3 ! noro 2448: gpmem_t av=avma;
1.1 noro 2449: GEN U;
2450: W=cgetg(lg(lo[i]),t_VECSMALL);
2451: for(k=1;k<lg(lo[i]);k++)
2452: W[k]=mael(lo,i,k)%f;
1.3 ! noro 2453: vecsmall_sort(W);
! 2454: U=vecsmall_uniq(W);
1.1 noro 2455: if (gegal(U, sg))
2456: {
2457: cgiv(U);
2458: V[j++]=lo[i];
2459: }
2460: else
2461: avma=av;
2462: }
2463: setlg(V,j);
2464: /*warning components of V point to W*/
2465: return gerepileupto(ltop,V);
2466: }
2467:
2468: static long
2469: galoisfrobeniustest(GEN aut, struct galois_lift *gl, GEN frob)
2470: {
1.3 ! noro 2471: gpmem_t ltop = avma;
1.1 noro 2472: GEN tlift = FpX_center(FpX_Fp_mul(aut,gl->den,gl->Q), gl->Q);
1.3 ! noro 2473: long res = poltopermtest(tlift, gl, frob);
! 2474: avma = ltop;
1.1 noro 2475: return res;
2476: }
2477:
2478: static GEN
2479: galoismakepsi(long g, GEN sg, GEN pf)
2480: {
2481: GEN psi=cgetg(g+1,t_VECSMALL);
2482: long i;
2483: for (i = 1; i < g; i++)
2484: psi[i] = sg[pf[i]];
2485: psi[g]=sg[1];
2486: return psi;
2487: }
2488:
2489: static GEN
1.3 ! noro 2490: galoisfrobeniuslift(GEN T, GEN den, GEN L, GEN Lden,
1.1 noro 2491: struct galois_frobenius *gf, struct galois_borne *gb)
2492: {
1.3 ! noro 2493: gpmem_t ltop=avma, av2;
1.1 noro 2494: struct galois_testlift gt;
2495: struct galois_lift gl;
2496: GEN res;
2497: long i,j,k;
2498: long n=lg(L)-1, deg=1, g=lg(gf->Tmod)-1;
2499: GEN F,Fp,Fe;
2500: GEN ip=stoi(gf->p), aut, frob;
2501: if (DEBUGLEVEL >= 4)
2502: fprintferr("GaloisConj:p=%ld deg=%ld fp=%ld\n", gf->p, deg, gf->fp);
2503: res = cgetg(lg(L), t_VECSMALL);
1.3 ! noro 2504: gf->psi = vecsmall_const(g,1);
1.1 noro 2505: av2=avma;
2506: initlift(T, den, ip, L, Lden, gb, &gl);
2507: aut = galoisdolift(&gl, res);
2508: if (!aut || galoisfrobeniustest(aut,&gl,res))
2509: {
2510: avma=av2;
2511: gf->deg = gf->fp;
2512: return res;
2513: }
2514: inittestlift(aut,gf->Tmod, &gl, >);
2515: gt.C=cgetg(gf->fp+1,t_VEC);
2516: for (i = 1; i <= gf->fp; i++)
2517: {
2518: gt.C[i]=lgetg(gt.g+1,t_VECSMALL);
2519: for(j = 1; j <= gt.g; j++)
2520: mael(gt.C,i,j)=0;
2521: }
2522: gt.Cd=gcopy(gt.C);
2523:
2524: F=factor(stoi(gf->fp));
2525: Fp=gtovecsmall((GEN)F[1]);
2526: Fe=gtovecsmall((GEN)F[2]);
2527: frob = cgetg(lg(L), t_VECSMALL);
2528: for(k=lg(Fp)-1;k>=1;k--)
2529: {
1.3 ! noro 2530: gpmem_t btop=avma;
1.1 noro 2531: GEN psi=NULL,fres=NULL,sg;
2532: long el=gf->fp, dg=1, dgf=1;
2533: long e,pr;
1.3 ! noro 2534: sg=perm_identity(1);
1.1 noro 2535: for(e=1;e<=Fe[k];e++)
2536: {
2537: long l;
2538: GEN lo;
2539: GEN pf;
2540: dg *= Fp[k]; el /= Fp[k];
2541: if ( DEBUGLEVEL>=4 )
2542: fprintferr("Trying degre %d.\n",dg);
2543: if (galoisfrobeniustest((GEN)gt.pauto[el+1],&gl,frob))
2544: {
2545: dgf = dg;
1.3 ! noro 2546: psi = vecsmall_const(g,1);
1.1 noro 2547: fres= gcopy(frob);
2548: continue;
2549: }
2550: disable_dbg(0);
2551: lo = listsousgroupes(dg, n / gf->fp);
2552: disable_dbg(-1);
2553: if (e!=1)
2554: lo = galoisfindgroups(lo, sg, dgf);
2555: if (DEBUGLEVEL >= 4)
2556: fprintferr("Galoisconj:Subgroups list:%Z\n", lo);
2557: for (l = 1; l < lg(lo); l++)
2558: if ( lg(lo[l])>2 &&
2559: frobeniusliftall((GEN)lo[l], el, &pf, &gl, >, frob))
2560: {
2561: sg = gcopy((GEN)lo[l]);
2562: psi = galoismakepsi(g,sg,pf);
2563: dgf = dg;
2564: fres=gcopy(frob);
2565: break;
2566: }
2567: if ( l == lg(lo) )
2568: break;
2569: }
2570: if (dgf==1) { avma=btop; continue; }
2571: pr=deg*dgf;
2572: if (deg==1)
2573: {
2574: for(i=1;i<lg(res);i++) res[i]=fres[i];
2575: for(i=1;i<lg(psi);i++) gf->psi[i]=psi[i];
2576: }
2577: else
2578: {
1.3 ! noro 2579: GEN cp=perm_mul(res,fres);
1.1 noro 2580: for(i=1;i<lg(res);i++) res[i]=cp[i];
2581: for(i=1;i<lg(psi);i++) gf->psi[i]=(dgf*gf->psi[i]+deg*psi[i])%pr;
2582: }
2583: deg=pr;
2584: avma=btop;
2585: }
2586: for (i = 1; i <= gf->fp; i++)
2587: for (j = 1; j <= gt.g; j++)
2588: if (mael(gt.C,i,j))
2589: gunclone(gmael(gt.C,i,j));
2590: if (DEBUGLEVEL>=4 && res)
2591: fprintferr("Best lift: %d\n",deg);
2592: if (deg==1)
2593: {
2594: avma=ltop;
2595: return NULL;
2596: }
2597: else
2598: {
2599: /*We need to normalise result so that psi[g]=1*/
2600: long im=itos(mpinvmod(stoi(gf->psi[g]),stoi(gf->fp)));
1.3 ! noro 2601: GEN cp=perm_pow(res, im);
1.1 noro 2602: for(i=1;i<lg(res);i++) res[i]=cp[i];
2603: for(i=1;i<lg(gf->psi);i++) gf->psi[i]=mulssmod(im,gf->psi[i],gf->fp);
2604: avma=av2;
2605: gf->deg=deg;
2606: return res;
2607: }
2608: }
2609:
2610: static GEN
1.3 ! noro 2611: galoisfindfrobenius(GEN T, GEN L, GEN den, struct galois_frobenius *gf,
1.1 noro 2612: struct galois_borne *gb, const struct galois_analysis *ga)
2613: {
1.3 ! noro 2614: gpmem_t lbot, ltop=avma;
! 2615: long Try=0;
1.1 noro 2616: long n = degpol(T), deg, gmask;
2617: byteptr primepointer = ga->primepointer;
2618: GEN Lden,frob;
2619: Lden=makeLden(L,den,gb);
2620: gf->deg=ga->deg;gf->p=ga->p; deg=ga->deg;
2621: gmask=(ga->group&ga_ext_2)?3:1;
2622: for (;;)
2623: {
1.3 ! noro 2624: gpmem_t av = avma;
1.1 noro 2625: long isram;
1.3 ! noro 2626: long i;
1.1 noro 2627: GEN ip,Tmod;
2628: ip = stoi(gf->p);
2629: Tmod = lift((GEN) factmod(T, ip));
2630: isram = 0;
2631: for (i = 1; i < lg(Tmod[2]) && !isram; i++)
2632: if (!gcmp1(gmael(Tmod,2,i)))
2633: isram = 1;
2634: if (isram == 0)
2635: {
2636: gf->fp = degpol(gmael(Tmod,1,1));
2637: for (i = 2; i < lg(Tmod[1]); i++)
2638: if (degpol(gmael(Tmod,1,i)) != gf->fp)
2639: {
2640: avma = ltop;
2641: return NULL; /* Not Galois polynomial */
2642: }
2643: lbot=avma;
2644: gf->Tmod=gcopy((GEN)Tmod[1]);
2645: if ( ((gmask&1) && gf->fp % deg == 0) || ((gmask&2) && gf->fp % 2== 0) )
2646: {
1.3 ! noro 2647: frob=galoisfrobeniuslift(T, den, L, Lden, gf, gb);
1.1 noro 2648: if (frob)
2649: {
2650: GEN *gptr[3];
2651: gptr[0]=&gf->Tmod;
2652: gptr[1]=&gf->psi;
2653: gptr[2]=&frob;
2654: gerepilemanysp(ltop,lbot,gptr,3);
2655: return frob;
2656: }
2657: if ((ga->group&ga_all_normal) && gf->fp % deg == 0)
2658: gmask&=~1;
2659: /*The first prime degree is always divisible by deg, so we don't
2660: * have to worry about ext_2 being used before regular supersolvable*/
2661: if (!gmask)
2662: {
2663: avma = ltop;
2664: return NULL;
2665: }
1.3 ! noro 2666: Try++;
! 2667: if ( (ga->group&ga_non_wss) && Try > n )
1.1 noro 2668: err(warner, "galoisconj _may_ hang up for this polynomial");
2669: }
2670: }
1.3 ! noro 2671: NEXT_PRIME_VIADIFF_CHECK(gf->p, primepointer);
1.1 noro 2672: if (DEBUGLEVEL >= 4)
2673: fprintferr("GaloisConj:next p=%ld\n", gf->p);
2674: avma = av;
2675: }
2676: }
2677:
2678: GEN
2679: galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne *gb,
1.3 ! noro 2680: const struct galois_analysis *ga);
! 2681: static GEN
! 2682: galoisgenfixedfield(GEN Tp, GEN Pmod, GEN V, GEN ip, struct galois_borne *gb, GEN Pg)
! 2683: {
! 2684: gpmem_t ltop=avma;
! 2685: GEN P, PL, Pden, PM, Pp, Pladicabs;
! 2686: GEN tau, PG;
! 2687: long g,gp;
! 2688: long x=varn(Tp);
! 2689: P=(GEN)V[2];
! 2690: PL=(GEN)V[1];
! 2691: gp=lg(Pmod)-1;
! 2692: Pp = FpX_red(P,ip);
! 2693: if (DEBUGLEVEL>=6)
! 2694: fprintferr("GaloisConj: Fixed field %Z\n",P);
! 2695: if (degpol(P)==2)
! 2696: {
! 2697: PG=cgetg(3,t_VEC);
! 2698: PG[1]=lgetg(2,t_VEC);
! 2699: PG[2]=lgetg(2,t_VECSMALL);
! 2700: mael(PG,1,1)=lgetg(3,t_VECSMALL);
! 2701: mael(PG,2,1)=2;
! 2702: mael3(PG,1,1,1)=2;
! 2703: mael3(PG,1,1,2)=1;
! 2704: tau=deg1pol(stoi(-1),negi((GEN)P[3]),x);
! 2705: tau = lift(gmul(tau,gmodulcp(gun,ip)));
! 2706: tau = FpX_FpXQ_compo((GEN) Pmod[gp], tau,Pp,ip);
! 2707: tau = FpX_gcd(Pp, tau,ip);
! 2708: tau = FpX_Fp_mul(tau,mpinvmod((GEN) tau[lgef(tau) - 1],ip),ip);
! 2709: for (g = 1; g <= gp; g++)
! 2710: if (gegal(tau, (GEN) Pmod[g]))
! 2711: break;
! 2712: if (g == lg(Pmod))
! 2713: return NULL;
! 2714: Pg[1]=g;
! 2715: }
! 2716: else
! 2717: {
! 2718: struct galois_analysis Pga;
! 2719: struct galois_borne Pgb;
! 2720: long j;
! 2721: galoisanalysis(P, &Pga, 0, 0);
! 2722: if (Pga.deg == 0)
! 2723: return NULL; /* Avoid computing the discriminant */
! 2724: Pgb.l = gb->l;
! 2725: Pden = galoisborne(P, NULL, &Pgb, Pga.ppp);
! 2726: Pladicabs=Pgb.ladicabs;
! 2727: if (Pgb.valabs > gb->valabs)
! 2728: {
! 2729: if (DEBUGLEVEL>=4)
! 2730: fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
! 2731: ,Pgb.valabs-gb->valabs);
! 2732: PL = rootpadicliftroots(P,PL,gb->l,Pgb.valabs);
! 2733: }
! 2734: PM = vandermondeinversemod(PL, P, Pden, Pgb.ladicabs);
! 2735: PG = galoisgen(P, PL, PM, Pden, &Pgb, &Pga);
! 2736: if (PG == gzero) return NULL;
! 2737: for (j = 1; j < lg(PG[1]); j++)
! 2738: {
! 2739: gpmem_t btop=avma;
! 2740: tau = permtopol(gmael(PG,1,j), PL, PM, Pden, Pladicabs, x);
! 2741: tau = lift(gmul(tau,gmodulcp(gun,ip)));
! 2742: tau = FpX_FpXQ_compo((GEN) Pmod[gp], tau,Pp,ip);
! 2743: tau = FpX_gcd(Pp, tau,ip);
! 2744: tau = FpX_Fp_mul(tau,mpinvmod((GEN) tau[lgef(tau) - 1],ip),ip);
! 2745: for (g = 1; g < lg(Pmod); g++)
! 2746: if (gegal(tau, (GEN) Pmod[g]))
! 2747: break;
! 2748: if (g == lg(Pmod))
! 2749: return NULL;
! 2750: avma=btop;
! 2751: Pg[j]=g;
! 2752: }
! 2753: }
! 2754: return gerepilecopy(ltop,PG);
! 2755: }
! 2756:
! 2757: GEN
! 2758: galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne *gb,
1.1 noro 2759: const struct galois_analysis *ga)
2760: {
2761: struct galois_test td;
2762: struct galois_frobenius gf;
1.3 ! noro 2763: gpmem_t lbot, ltop2, ltop = avma;
1.1 noro 2764: long n, p, deg;
1.3 ! noro 2765: long fp;
1.1 noro 2766: long x;
2767: int i, j;
2768: GEN Lden, sigma;
2769: GEN Tmod, res, pf = gzero, split, ip;
2770: GEN frob;
2771: GEN O;
1.3 ! noro 2772: GEN PG, Pg;
1.1 noro 2773: n = degpol(T);
2774: if (!ga->deg)
2775: return gzero;
2776: x = varn(T);
2777: if (DEBUGLEVEL >= 9)
2778: fprintferr("GaloisConj:denominator:%Z\n", den);
2779: if (n == 12 && ga->ord==3) /* A4 is very probable,so test it first */
2780: {
1.3 ! noro 2781: gpmem_t av = avma;
1.1 noro 2782: if (DEBUGLEVEL >= 4)
2783: fprintferr("GaloisConj:Testing A4 first\n");
2784: inittest(L, M, gb->bornesol, gb->ladicsol, &td);
2785: lbot = avma;
2786: PG = a4galoisgen(T, &td);
2787: freetest(&td);
2788: if (PG != gzero)
2789: return gerepile(ltop, lbot, PG);
2790: avma = av;
2791: }
2792: if (n == 24 && ga->ord==3) /* S4 is very probable,so test it first */
2793: {
1.3 ! noro 2794: gpmem_t av = avma;
1.1 noro 2795: struct galois_lift gl;
2796: if (DEBUGLEVEL >= 4)
2797: fprintferr("GaloisConj:Testing S4 first\n");
2798: lbot = avma;
2799: Lden=makeLden(L,den,gb);
2800: initlift(T, den, stoi(ga->p4), L, Lden, gb, &gl);
2801: PG = s4galoisgen(&gl);
2802: if (PG != gzero)
2803: return gerepile(ltop, lbot, PG);
2804: avma = av;
2805: }
1.3 ! noro 2806: frob=galoisfindfrobenius(T, L, den, &gf, gb, ga);
1.1 noro 2807: if (!frob)
2808: {
2809: ltop=avma;
2810: return gzero;
2811: }
1.3 ! noro 2812: p=gf.p;
1.1 noro 2813: ip=stoi(p);
2814: Tmod=gf.Tmod;
1.3 ! noro 2815: fp=gf.fp;
! 2816: O = perm_cycles(frob);
1.1 noro 2817: split = splitorbite(O);
2818: deg=lg(O[1])-1;
2819: sigma = permtopol(frob, L, M, den, gb->ladicabs, x);
2820: if (DEBUGLEVEL >= 9)
2821: fprintferr("GaloisConj:Orbite:%Z\n", O);
2822: if (deg == n) /* Cyclique */
2823: {
2824: lbot = avma;
2825: res = cgetg(3, t_VEC);
2826: res[1] = lgetg(lg(split[1]), t_VEC);
2827: res[2] = lgetg(lg(split[2]), t_VECSMALL);
2828: for (i = 1; i < lg(split[1]); i++)
2829: {
2830: mael(res,1,i) = lcopy(gmael(split,1,i));
2831: mael(res,2,i) = mael(split,2,i);
2832: }
2833: return gerepile(ltop, lbot, res);
2834: }
2835: if (DEBUGLEVEL >= 9)
2836: fprintferr("GaloisConj:Frobenius:%Z\n", sigma);
1.3 ! noro 2837: Pg=cgetg(lg(O),t_VECSMALL);
1.1 noro 2838: {
1.3 ! noro 2839: gpmem_t btop=avma;
! 2840: GEN V, sym, dg, Tp, Sp, Pmod;
1.1 noro 2841: sym=cgetg(lg(L),t_VECSMALL);
2842: dg=cgetg(lg(L),t_VECSMALL);
2843: V = fixedfieldsympol(O, L, gb->ladicabs, gb->l, ip, sym, dg, x);
2844: Tp = FpX_red(T,ip);
2845: Sp = fixedfieldpolsigma(sigma,ip,Tp,sym,dg,deg);
2846: Pmod = fixedfieldfactmod(Sp,ip,Tmod);
1.3 ! noro 2847: PG=galoisgenfixedfield(Tp, Pmod, V, ip, gb, Pg);
! 2848: if (PG == NULL)
1.1 noro 2849: {
1.3 ! noro 2850: avma = ltop;
! 2851: return gzero;
1.1 noro 2852: }
1.3 ! noro 2853: if (DEBUGLEVEL >= 4)
! 2854: fprintferr("GaloisConj:Back to Earth:%Z\n", PG);
! 2855: PG=gerepileupto(btop, PG);
1.1 noro 2856: }
2857: inittest(L, M, gb->bornesol, gb->ladicsol, &td);
2858: lbot = avma;
2859: res = cgetg(3, t_VEC);
2860: res[1] = lgetg(lg(PG[1]) + lg(split[1]) - 1, t_VEC);
2861: res[2] = lgetg(lg(PG[1]) + lg(split[1]) - 1, t_VECSMALL);
2862: for (i = 1; i < lg(split[1]); i++)
2863: {
2864: ((GEN **) res)[1][i] = gcopy(((GEN **) split)[1][i]);
2865: ((GEN *) res)[2][i] = ((GEN *) split)[2][i];
2866: }
2867: for (i = lg(split[1]); i < lg(res[1]); i++)
2868: ((GEN **) res)[1][i] = cgetg(n + 1, t_VECSMALL);
2869: ltop2 = avma;
2870: for (j = 1; j < lg(PG[1]); j++)
2871: {
1.3 ! noro 2872: GEN B;
! 2873: long t;
1.1 noro 2874: long w, sr, dss;
2875: if (DEBUGLEVEL >= 6)
2876: fprintferr("GaloisConj:G[%d]=%Z d'ordre relatif %d\n", j,
2877: ((GEN **) PG)[1][j], ((long **) PG)[2][j]);
1.3 ! noro 2878: B = perm_cycles(gmael(PG,1,j));
1.1 noro 2879: if (DEBUGLEVEL >= 6)
2880: fprintferr("GaloisConj:B=%Z\n", B);
1.3 ! noro 2881: w = gf.psi[Pg[j]];
1.1 noro 2882: dss = deg / cgcd(deg, w - 1);
2883: sr = 1;
2884: for (i = 1; i < lg(B[1]) - 1; i++)
2885: sr = (1 + sr * w) % deg;
2886: sr = cgcd(sr, deg);
2887: if (DEBUGLEVEL >= 6)
2888: fprintferr("GaloisConj:w=%ld [%ld] sr=%ld dss=%ld\n", w, deg, sr, dss);
2889: for (t = 0; t < sr; t += dss)
2890: {
2891: pf = testpermutation(O, B, w, t, sr, &td);
2892: if (pf != gzero)
2893: break;
2894: }
2895: if (pf == gzero)
2896: {
2897: freetest(&td);
2898: avma = ltop;
2899: return gzero;
2900: }
2901: else
2902: {
2903: int i;
2904: for (i = 1; i <= n; i++)
2905: ((GEN **) res)[1][lg(split[1]) + j - 1][i] = pf[i];
2906: ((GEN *) res)[2][lg(split[1]) + j - 1] = ((GEN *) PG)[2][j];
2907: }
2908: avma = ltop2;
2909: }
2910: if (DEBUGLEVEL >= 4)
2911: fprintferr("GaloisConj:Fini!\n");
2912: freetest(&td);
2913: return gerepile(ltop, lbot, res);
2914: }
2915:
2916: /*
2917: * T: polynome ou nf den optionnel: si présent doit etre un multiple du
2918: * dénominateur commun des solutions Si T est un nf et den n'est pas présent,
2919: * le denominateur de la base d'entiers est pris pour den
2920: */
2921: GEN
2922: galoisconj4(GEN T, GEN den, long flag, long karma)
2923: {
1.3 ! noro 2924: gpmem_t ltop = avma;
1.1 noro 2925: GEN G, L, M, res, aut, grp=NULL;/*keep gcc happy on the wall*/
2926: struct galois_analysis ga;
2927: struct galois_borne gb;
2928: int n, i, j, k;
2929: if (typ(T) != t_POL)
2930: {
2931: T = checknf(T);
1.3 ! noro 2932: if (!den) den = Q_denom((GEN)T[7]);
1.1 noro 2933: T = (GEN) T[1];
2934: }
2935: n = degpol(T);
2936: if (n <= 0)
2937: err(constpoler, "galoisconj4");
2938: for (k = 2; k <= n + 2; k++)
2939: if (typ(T[k]) != t_INT)
2940: err(talker, "polynomial not in Z[X] in galoisconj4");
2941: if (!gcmp1((GEN) T[n + 2]))
2942: err(talker, "non-monic polynomial in galoisconj4");
2943: n = degpol(T);
2944: if (n == 1) /* Too easy! */
2945: {
2946: if (!flag)
2947: {
2948: res = cgetg(2, t_COL);
2949: res[1] = (long) polx[varn(T)];
2950: return res;
2951: }
2952: ga.l = 3;
2953: ga.deg = 1;
2954: ga.ppp = 1;
2955: den = gun;
2956: }
2957: else
2958: galoisanalysis(T, &ga, 1, karma);
2959: if (ga.deg == 0)
2960: {
2961: avma = ltop;
2962: return stoi(ga.p); /* Avoid computing the discriminant */
2963: }
2964: if (den)
2965: {
2966: if (typ(den) != t_INT)
2967: err(talker, "Second arg. must be integer in galoisconj4");
2968: den = absi(den);
2969: }
2970: gb.l = stoi(ga.l);
1.3 ! noro 2971: if (DEBUGLEVEL >= 1) (void)timer2();
1.1 noro 2972: den = galoisborne(T, den, &gb, ga.ppp);
2973: if (DEBUGLEVEL >= 1)
2974: msgtimer("galoisborne()");
2975: L = rootpadicfast(T, gb.l, gb.valabs);
2976: if (DEBUGLEVEL >= 1)
2977: msgtimer("rootpadicfast()");
2978: M = vandermondeinversemod(L, T, den, gb.ladicabs);
2979: if (DEBUGLEVEL >= 1)
2980: msgtimer("vandermondeinversemod()");
2981: if (n == 1)
2982: {
2983: G = cgetg(3, t_VEC);
2984: G[1] = lgetg(1, t_VECSMALL);
2985: G[2] = lgetg(1, t_VECSMALL);
2986: }
2987: else
2988: G = galoisgen(T, L, M, den, &gb, &ga);
2989: if (DEBUGLEVEL >= 6)
2990: fprintferr("GaloisConj:%Z\n", G);
2991: if (G == gzero)
2992: {
2993: avma = ltop;
2994: return gzero;
2995: }
1.3 ! noro 2996: if (DEBUGLEVEL >= 1) (void)timer2();
1.1 noro 2997: if (flag)
2998: {
2999: grp = cgetg(9, t_VEC);
3000: grp[1] = (long) gcopy(T);
3001: grp[2] = lgetg(4,t_VEC); /*Make K.B. happy(8 components)*/
3002: mael(grp,2,1) = lstoi(ga.l);
3003: mael(grp,2,2) = lstoi(gb.valabs);
3004: mael(grp,2,3) = licopy(gb.ladicabs);
3005: grp[3] = (long) gcopy(L);
3006: grp[4] = (long) gcopy(M);
3007: grp[5] = (long) gcopy(den);
3008: grp[7] = (long) gcopy((GEN) G[1]);
3009: grp[8] = (long) gcopy((GEN) G[2]);
3010: }
3011: res = cgetg(n + 1, t_VEC);
1.3 ! noro 3012: res[1] = (long) perm_identity(n);
1.1 noro 3013: k = 1;
3014: for (i = 1; i < lg(G[1]); i++)
3015: {
3016: int c = k * (((long **) G)[2][i] - 1);
3017: for (j = 1; j <= c; j++) /* I like it */
1.3 ! noro 3018: res[++k] = (long) perm_mul((GEN) res[j], ((GEN **) G)[1][i]);
1.1 noro 3019: }
3020: if (flag)
3021: {
3022: grp[6] = (long) res;
3023: return gerepileupto(ltop, grp);
3024: }
3025: aut = galoisgrouptopol(res,L,M,den,gb.ladicsol, varn(T));
3026: if (DEBUGLEVEL >= 1)
3027: msgtimer("Calcul polynomes");
1.3 ! noro 3028: return gerepileupto(ltop, gen_sort(aut, 0, cmp_pol));
1.1 noro 3029: }
3030:
3031: /* Calcule le nombre d'automorphisme de T avec forte probabilité */
3032: /* pdepart premier nombre premier a tester */
3033: long
3034: numberofconjugates(GEN T, long pdepart)
3035: {
1.3 ! noro 3036: gpmem_t ltop2, ltop = avma;
1.1 noro 3037: long n, p, nbmax, nbtest;
3038: long card;
3039: byteptr primepointer;
3040: int i;
3041: GEN L;
3042: n = degpol(T);
3043: card = sturm(T);
3044: card = cgcd(card, n - card);
3045: nbmax = (n<<1) + 1;
3046: if (nbmax < 20) nbmax=20;
3047: nbtest = 0;
3048: L = cgetg(n + 1, t_VECSMALL);
3049: ltop2 = avma;
3050: for (p = 0, primepointer = diffptr; nbtest < nbmax && card > 1;)
3051: {
1.3 ! noro 3052: long s;
1.1 noro 3053: long isram;
3054: GEN S;
1.3 ! noro 3055:
! 3056: NEXT_PRIME_VIADIFF_CHECK(p,primepointer);
1.1 noro 3057: if (p < pdepart)
3058: continue;
3059: S = (GEN) simplefactmod(T, stoi(p));
3060: isram = 0;
3061: for (i = 1; i < lg(S[2]) && !isram; i++)
3062: if (!gcmp1(((GEN **) S)[2][i]))
3063: isram = 1;
3064: if (isram == 0)
3065: {
3066: for (i = 1; i <= n; i++)
3067: L[i] = 0;
3068: for (i = 1; i < lg(S[1]) && !isram; i++)
3069: L[itos(((GEN **) S)[1][i])]++;
3070: s = L[1];
3071: for (i = 2; i <= n; i++)
3072: s = cgcd(s, L[i] * i);
3073: card = cgcd(s, card);
3074: }
3075: if (DEBUGLEVEL >= 6)
3076: fprintferr("NumberOfConjugates:Nbtest=%ld,card=%ld,p=%ld\n", nbtest,
3077: card, p);
3078: nbtest++;
3079: avma = ltop2;
3080: }
3081: if (DEBUGLEVEL >= 2)
3082: fprintferr("NumberOfConjugates:card=%ld,p=%ld\n", card, p);
3083: avma = ltop;
3084: return card;
3085: }
3086:
3087: GEN
3088: galoisconj0(GEN nf, long flag, GEN d, long prec)
3089: {
1.3 ! noro 3090: gpmem_t ltop;
1.1 noro 3091: GEN G, T;
3092: long card;
3093: if (typ(nf) != t_POL)
3094: {
3095: nf = checknf(nf);
3096: T = (GEN) nf[1];
3097: }
3098: else
3099: T = nf;
3100: switch (flag)
3101: {
3102: case 0:
3103: ltop = avma;
3104: G = galoisconj4(nf, d, 0, 0);
3105: if (typ(G) != t_INT) /* Success */
3106: return G;
3107: else
3108: {
3109: card = numberofconjugates(T, G == gzero ? 2 : itos(G));
3110: avma = ltop;
3111: if (card != 1)
3112: {
3113: if (typ(nf) == t_POL)
3114: {
3115: G = galoisconj2pol(nf, card, prec);
3116: if (lg(G) <= card)
3117: err(warner, "conjugates list may be incomplete in nfgaloisconj");
3118: return G;
3119: }
3120: else
3121: return galoisconj(nf);
3122: }
3123: }
3124: break; /* Failure */
3125: case 1:
3126: return galoisconj(nf);
3127: case 2:
3128: return galoisconj2(nf, degpol(T), prec);
3129: case 4:
3130: G = galoisconj4(nf, d, 0, 0);
3131: if (typ(G) != t_INT)
3132: return G;
3133: break; /* Failure */
3134: default:
3135: err(flagerr, "nfgaloisconj");
3136: }
3137: G = cgetg(2, t_COL);
3138: G[1] = (long) polx[varn(T)];
3139: return G; /* Failure */
3140: }
3141:
3142:
3143:
3144: /******************************************************************************/
3145: /* Isomorphism between number fields */
3146: /******************************************************************************/
3147: long
3148: isomborne(GEN P, GEN den, GEN p)
3149: {
1.3 ! noro 3150: gpmem_t ltop=avma;
1.1 noro 3151: struct galois_borne gb;
3152: gb.l=p;
1.3 ! noro 3153: (void)galoisborne(P,den,&gb,degpol(P));
1.1 noro 3154: avma=ltop;
3155: return gb.valsol;
3156: }
3157:
3158:
3159:
3160: /******************************************************************************/
3161: /* Galois theory related algorithms */
3162: /******************************************************************************/
3163: GEN
3164: checkgal(GEN gal)
3165: {
3166: if (typ(gal) == t_POL)
3167: err(talker, "please apply galoisinit first");
3168: if (typ(gal) != t_VEC || lg(gal) != 9)
3169: err(talker, "Not a Galois field in a Galois related function");
3170: return gal;
3171: }
3172:
3173: GEN
3174: galoisinit(GEN nf, GEN den, long karma)
3175: {
3176: GEN G;
3177: G = galoisconj4(nf, den, 1, karma);
3178: if (typ(G) == t_INT)
3179: err(talker,
3180: "galoisinit: field not Galois or Galois group not weakly super solvable");
3181: return G;
3182: }
3183:
3184: GEN
3185: galoispermtopol(GEN gal, GEN perm)
3186: {
3187: GEN v;
3188: long t = typ(perm);
3189: int i;
3190: gal = checkgal(gal);
3191: switch (t)
3192: {
3193: case t_VECSMALL:
3194: return permtopol(perm, (GEN) gal[3], (GEN) gal[4], (GEN) gal[5],
3195: gmael(gal,2,3), varn((GEN) gal[1]));
3196: case t_VEC:
3197: case t_COL:
3198: case t_MAT:
3199: v = cgetg(lg(perm), t);
3200: for (i = 1; i < lg(v); i++)
3201: v[i] = (long) galoispermtopol(gal, (GEN) perm[i]);
3202: return v;
3203: }
3204: err(typeer, "galoispermtopol");
3205: return NULL; /* not reached */
3206: }
3207:
3208:
3209: GEN
3210: galoiscosets(GEN O, GEN perm)
3211: {
3212: long i,j,k,u;
3213: long o = lg(O)-1, f = lg(O[1])-1;
3214: GEN C = cgetg(lg(O),t_VECSMALL);
1.3 ! noro 3215: gpmem_t av=avma;
1.1 noro 3216: GEN RC=cgetg(lg(perm),t_VECSMALL);
3217: for(i=1;i<lg(RC);i++)
3218: RC[i]=0;
3219: u=mael(O,1,1);
3220: for(i=1,j=1;j<=o;i++)
3221: {
3222: if (RC[mael(perm,i,u)])
3223: continue;
3224: for(k=1;k<=f;k++)
3225: RC[mael(perm,i,mael(O,1,k))]=1;
3226: C[j++]=i;
3227: }
3228: avma=av;
3229: return C;
3230: }
3231:
3232: GEN
3233: fixedfieldfactor(GEN L, GEN O, GEN perm, GEN M, GEN den, GEN mod,
3234: long x,long y)
3235: {
1.3 ! noro 3236: gpmem_t ltop=avma;
1.1 noro 3237: GEN F,G,V,res,cosets;
3238: int i, j, k;
3239: F=cgetg(lg(O[1])+1,t_COL);
3240: F[lg(O[1])]=un;
3241: G=cgetg(lg(O),t_VEC);
3242: for (i = 1; i < lg(O); i++)
3243: {
3244: GEN Li = cgetg(lg(O[i]),t_VEC);
3245: for (j = 1; j < lg(O[i]); j++)
3246: Li[j] = L[mael(O,i,j)];
3247: G[i]=(long)FpV_roots_to_pol(Li,mod,x);
3248: }
3249:
3250: cosets=galoiscosets(O,perm);
3251: if (DEBUGLEVEL>=4) fprintferr("GaloisFixedField:cosets=%Z \n",cosets);
3252: V=cgetg(lg(O),t_COL);
3253: if (DEBUGLEVEL>=6) fprintferr("GaloisFixedField:den=%Z mod=%Z \n",den,mod);
3254: res=cgetg(lg(O),t_VEC);
3255: for (i = 1; i < lg(O); i++)
3256: {
1.3 ! noro 3257: gpmem_t av=avma;
1.1 noro 3258: long ii,jj;
3259: GEN G=cgetg(lg(O),t_VEC);
3260: for (ii = 1; ii < lg(O); ii++)
3261: {
3262: GEN Li = cgetg(lg(O[ii]),t_VEC);
3263: for (jj = 1; jj < lg(O[ii]); jj++)
3264: Li[jj] = L[mael(perm,cosets[i],mael(O,ii,jj))];
3265: G[ii]=(long)FpV_roots_to_pol(Li,mod,x);
3266: }
3267: for (j = 1; j < lg(O[1]); j++)
3268: {
3269: for(k = 1; k < lg(O); k++)
3270: V[k]=mael(G,k,j+1);
3271: F[j]=(long) vectopol(V, M, den, mod, y);
3272: }
3273: res[i]=(long)gerepileupto(av,gtopolyrev(F,x));
3274: }
3275: return gerepileupto(ltop,res);
3276: }
3277:
3278: GEN
3279: galoisfixedfield(GEN gal, GEN perm, long flag, long y)
3280: {
1.3 ! noro 3281: gpmem_t lbot, ltop = avma;
1.1 noro 3282: GEN L, P, S, PL, O, res, mod;
3283: long x, n;
3284: int i;
3285: gal = checkgal(gal);
3286: x = varn((GEN) gal[1]);
3287: L = (GEN) gal[3]; n=lg(L)-1;
3288: mod = gmael(gal,2,3);
3289: if (flag<0 || flag>2)
3290: err(flagerr, "galoisfixedfield");
3291: if (typ(perm) == t_VEC)
3292: {
1.3 ! noro 3293: for (i = 1; i < lg(perm); i++)
! 3294: if (typ(perm[i]) != t_VECSMALL || lg(perm[i])!=n+1)
! 3295: err(typeer, "galoisfixedfield");
! 3296: O = vecperm_orbits(perm, n);
1.1 noro 3297: }
3298: else if (typ(perm) != t_VECSMALL || lg(perm)!=n+1 )
1.3 ! noro 3299: {
1.1 noro 3300: err(typeer, "galoisfixedfield");
1.3 ! noro 3301: return NULL; /* not reached */
! 3302: }
! 3303: else
! 3304: O = perm_cycles(perm);
1.1 noro 3305: {
3306: GEN sym=cgetg(lg(L),t_VECSMALL);
3307: GEN dg=cgetg(lg(L),t_VECSMALL);
3308: GEN V;
3309: V = fixedfieldsympol(O, L, mod, gmael(gal,2,1), gun, sym, dg, x);
3310: P=(GEN)V[2];
3311: PL=(GEN)V[1];
3312: }
3313: if (flag==1)
3314: return gerepileupto(ltop,P);
3315: S = fixedfieldinclusion(O, PL);
3316: S = vectopol(S, (GEN) gal[4], (GEN) gal[5], mod, x);
3317: if (flag==0)
3318: {
3319: lbot = avma;
3320: res = cgetg(3, t_VEC);
3321: res[1] = (long) gcopy(P);
3322: res[2] = (long) gmodulcp(S, (GEN) gal[1]);
3323: return gerepile(ltop, lbot, res);
3324: }
3325: else
3326: {
3327: GEN PM,Pden;
3328: {
3329: struct galois_borne Pgb;
3330: long val=itos(gmael(gal,2,2));
3331: Pgb.l = gmael(gal,2,1);
3332: Pden = galoisborne(P, NULL, &Pgb, 1);
3333: if (Pgb.valabs > val)
3334: {
1.3 ! noro 3335: if (DEBUGLEVEL>=4)
! 3336: fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
! 3337: ,Pgb.valabs-val);
! 3338: PL = rootpadicliftroots(P,PL,Pgb.l,Pgb.valabs);
! 3339: L = rootpadicliftroots((GEN) gal[1],L,Pgb.l,Pgb.valabs);
! 3340: mod = Pgb.ladicabs;
1.1 noro 3341: }
3342: }
3343: PM = vandermondeinversemod(PL, P, Pden, mod);
3344: lbot = avma;
3345: if (y==-1)
3346: y = fetch_user_var("y");
3347: if (y<=x)
3348: err(talker,"priority of optional variable too high in galoisfixedfield");
3349: res = cgetg(4, t_VEC);
3350: res[1] = (long) gcopy(P);
3351: res[2] = (long) gmodulcp(S, (GEN) gal[1]);
3352: res[3] = (long) fixedfieldfactor(L,O,(GEN)gal[6],
1.3 ! noro 3353: PM,Pden,mod,x,y);
1.1 noro 3354: return gerepile(ltop, lbot, res);
3355: }
3356: }
3357:
3358: /*return 1 if gal is abelian, else 0*/
3359: long
3360: galoistestabelian(GEN gal)
3361: {
1.3 ! noro 3362: gpmem_t ltop=avma;
! 3363: GEN G;
! 3364: long test;
1.1 noro 3365: gal = checkgal(gal);
1.3 ! noro 3366: G=cgetg(3,t_VEC);
! 3367: G[1]=gal[7];
! 3368: G[2]=gal[8];
! 3369: test=group_isabelian(G);
! 3370: avma=ltop;
! 3371: return test;
1.1 noro 3372: }
3373:
3374: GEN galoisisabelian(GEN gal, long flag)
3375: {
3376: long i, j;
1.3 ! noro 3377: long test, n;
1.1 noro 3378: GEN M;
3379: gal = checkgal(gal);
3380: test=galoistestabelian(gal);
3381: if (!test) return gzero;
3382: if (flag==1) return gun;
3383: if (flag) err(flagerr,"galoisisabelian");
1.3 ! noro 3384: n=lg(gal[7]);
1.1 noro 3385: M=cgetg(n,t_MAT);
3386: for(i=1;i<n;i++)
3387: {
1.3 ! noro 3388: gpmem_t btop;
1.1 noro 3389: GEN P;
3390: long k;
3391: M[i]=lgetg(n,t_COL);
3392: btop=avma;
1.3 ! noro 3393: P=perm_pow(gmael(gal,7,i),mael(gal,8,i));
1.1 noro 3394: for(j=1;j<lg(gal[6]);j++)
3395: if (gegal(P,gmael(gal,6,j)))
3396: break;
3397: avma=btop;
3398: if (j==lg(gal[6])) err(talker,"wrong argument in galoisisabelian");
3399: j--;
3400: for(k=1;k<i;k++)
3401: {
3402: mael(M,i,k)=lstoi(j%mael(gal,8,k));
3403: j/=mael(gal,8,k);
3404: }
3405: mael(M,i,k++)=lstoi(mael(gal,8,i));
3406: for( ;k<n;k++)
3407: mael(M,i,k)=zero;
3408: }
3409: return M;
3410: }
1.3 ! noro 3411: GEN galoissubgroups(GEN G)
1.1 noro 3412: {
1.3 ! noro 3413: gpmem_t ltop=avma;
! 3414: GEN H;
! 3415: G = checkgal(G);
! 3416: H=cgetg(3,t_VEC);
! 3417: H[1]=G[7];
! 3418: H[2]=G[8];
! 3419: return gerepileupto(ltop,group_subgroups(H));
1.1 noro 3420: }
3421:
1.3 ! noro 3422: GEN galoissubfields(GEN G, long flag, long v)
1.1 noro 3423: {
1.3 ! noro 3424: gpmem_t ltop=avma;
1.1 noro 3425: long i;
1.3 ! noro 3426: GEN L = galoissubgroups(G);
! 3427: long l2 = lg(L);
! 3428: GEN p3 = cgetg(l2, t_VEC);
! 3429: for (i = 1; i < l2; ++i)
! 3430: p3[i] = (long) galoisfixedfield(G, gmael(L,i,1), flag, v);
! 3431: return gerepileupto(ltop,p3);
1.1 noro 3432: }
3433:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to galconj.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath