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